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and  a quadratic  constraint  imposed  on  the  variations  of  the 
control,  the  parameter,  and  the  missing  components  of  the 
initial  state.  The  restoration  phase  involves  one  or  more 
iterations:  in  each  iteration,  the  norm  squared  of  the 
variations  of  the  control,  the  parameter,  and  the  missing 
components  of  the  initial  state  is  minimized,  subject  to 
the  linearized  constraints. 

The  algorithm  has  two  principal  properties:  (i)  it 
provides  a sequence  of  suboptimal  solutions;  the  functions 
obtained  at  the  end  of  each  cycle  satisfy  the  constraints  to 
a predetermined  accuracy;  and  (ii)  the  value  of  the 
functional  at  the  end  of  any  complete  conjugate  gradient- 
restoration  cycle  is  smaller  than  the  value  of  the  same 
functional  at  the  beginning  of  that  cycle. 

Th’e  sequential  conjugate  gradient-restoration  algorithm 
presented  hero  differs  from  previous  work,  in  that  it  is  not 
required  that  the  state  vector  be  given  at  the  initial  point. 
Instead,  the  initial  conditions  can  be  absolutely  general. 
Since  the  present  algorithm  is  capable  of  handling  general 
final  conditions,  it  is  suitable  for  the  solution  of  optimal 
control  problems  with  general  boundary  conditions.  Its 
importance  lies  in  the  fact  that  many  optimal  control  prob- 
lems can  bd  reduced  to  the  format  considered  here.  This 
includes  (i)  problems  with  state  equality  constraints,  (ii) 
problems  with  control  equality  constraints,  (iii)  problems 
with  state-derivative  equality  constraints,  (iv)  problems 
with  state  inequality  constraints,  (v)  problems  with  control 
inequality  constraints,  (vi)  problems  with  state-derivative 
inequality  constraints,  and  (vii)  minimax  problems  of  optimal 
control . 

7 Twelve  numerical  examples  are  presented  in  order  to 
illustrate  the  performance  of  the  algorithm.  The  numerical 
results  show  the  feasibility  as  well  as  the  convergence 
characteristics  of  the  present  algorithm. 
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for  Optimal  Control  Problems 
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Abstract.  This  paper  considers  the  numerical  solution  o 


the  problem  of  minimizing  a functional  I subject  to  difforen 


tial  constraints,  nondifferential  constraints 


boundary  conditions.  It  consists  of  finding  the  state  x(t) 


the  control  u(t),  and  the  parameter  u 


o that  the  functional 


I is  minimized  while  the  constraints  and  the  boundary  condi 
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The  approach  taken  is  a sequence  of  two-phase  cycles, 
composed  of  a conjugate  gradient  phase  and  a restoration 
phase.  The  conjugate  gradient  phase  involves  one  iteration 
and  is  designed  to  decrease  the  value  of  the  functional, 
while  the  constraints  are  satisfied  to  first  order.  During 
this  iteration,  the  first  variation  of  the  functional  is 
minimized,  subject  to  linearized  constraints.  The  minimiza- 
tion is  performed  over  the  class  of  variations  of  the  control, 
the  parameter,  and  the  missing  components  of  the  initial 
state  which  are  equidistant  from  some  constant  multiple  of 
the  corresponding  variations  of  the  previous  conjugate  gra- 
dient phase.  The  sequence  of  conjugate  gradient  phases 
generated  by  the  algorithm  is  such  that,  for  the  special  case  of  a 
quadratic  functional  subject  to  linear  constraints,  various 
orthogonality  and  conjugacy  conditions  hold.  The  restoration 
phase  involves  one  or  more  iterations  and  is  designed  to 
force  constraint  satisfaction  to  a predetermined  accuracy, 
while  the  norm  squared  of  the  variations  of  the  control,  the 
parameter,  and  the  missing  components  of  the  initial  state 
is  minimized. 

The  principal  property  of  the  algorithm  is  that  it  pro- 
duces a sequence  of  feasible  suboptimal  solutions:  the  fun- 
ctions obtained  at  the  end  of  each  cycle  satisfy  the  cons- 
traints to  a predetermined  accuracy.  Therefore,  the  values 


sequence  are  comparable. 

The  stepsize  of  the  conjugate  gradient  phase  is  deter- 
mined by  a one-dimensional  search  on  the  augmented  functional 
J,  while  the  stepsize  of  the  restoration  phase  is  obtained 
by  a one-dimensional  search  on  the  constraint  error  P.  The 


conjugate  gradient  stepsize  and  the  restoration  stepsize  are 
chosen  so  that  the  restoration  phase  preserves  the  descent 


property  of  the  conjugate  gradient  phase.  Therefore,  the 


value  of  the  functional  I at  the  end  of  any  complete  conju- 
gate gradient-restoration  cycle  is  smaller  than  the  value  of 
the  same  functional  at  the  beginning  of  that  cycle.  Of 
course,  restarting  the  algorithm  might  be  occasionally  neces- 
sary . 

The  sequential  conjugate  gradient-restoration  algorithm 
presented  here  differs  from  that  of  Refs,  3 and  4,  in  that  it 
is  not  required  that  the  state  vector  be  given  at  the  initial 
point.  Instead,  the  initial  conditions  can  be  absolutely 
general.  In  analogy  with  Refs.  3 and  4,  the  present  algorithm 
is  capable  of  handling  general  final  conditions;  therefore,  it 
is  suitable  for  the  solution  of  optimal  control  problems 
with  general  boundary  conditions.  Its  importance  lies  in  the 
fact  that  many  optimal  control  problems  involve  initial  con- 
ditions of  the  type  considered  here. 
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Twelve  numerical  examples  are  presented  in  order  to 
illustrate  the  performance  of  the  algorithm.  The  numerical 
results  show  the  feasibility  as  well  as  the  convergence 
characteristics  of  the  present  algorithm. 

Key  Words.  Optimal  control,  numerical  methods,  computing 
methods,  gradient  methods,  gradient-restoration  algorithms, 
sequential  gradient-restoration  algorithms,  conjugate  gradient- 
restoration  algorithms,  sequential  conjugate  gradient- 
restoration  algorithms,  nondifferential  constraints,  bounded 
control,  bounded  state,  general  boundary  conditions. 
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1 . Introduction 

Over  the  past  several  years,  a successful  family  of  al- 
gorithms for  the  solution  of  optimal  control  problems  involv- 
ing differential  constraints,  nondifferential  constraints, 
and  terminal  constraints  has  been  developed  at  Rice  Universi- 
ty by  Miele  and  his  associates  (see  Refs.  1-4).  They  are 
known  as  sequential  gradient-restoration  algorithms  and  have 
been  designed  for  the  solution  of  different  classes  of  optimal 
control  problems.  Some  of  these  algorithms  are  of  the 
ordinary-gradient  type  (Refs.  1-2),  while  the  rest  are  of  the 
conjugate-gradient  type  (Refs.  2-4).  All  of  them  have  shown 
to  be  robust  and  reliable;  however,  they  all  require  the 
state  vector  to  be  given  at  the  initial  point.  Owing  to  the 
fact  that  optimal  control  problems  exist,  which  require  satis- 
faction of  more  general  boundary  conditions,  the  task  of  ex- 
tending the  aforementioned  family  of  algorithms  must  be 
undertaken,  with  these  ideas  in  mind:  (i)  to  retain  the 
robustness,  reliability,  and  convergence  characteristics  of 
the  algorithms  discussed  in  Refs.  1-4;  (ii)  to  be  able  to 
handle  all  of  the  optimal  control  problems  treated  in  Refs. 
1-4;  and  (iii)  to  have  the  additional  capability  of  handling 
optimal  control  problems  with  general  boundary  conditions. 

In  Ref.  5,  an  algorithm  of  the  ordinary  gradient  type 
was  developed,  extending  the  work  of  Ref.  1 to  problems  with 
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general  boundary  conditions.  In  this  report,  an  algorithm  of 
the  conjugate  gradient  type  is  developed,  extending  the  work 
of  Refs.  3-4  to  problems  with  general  boundary  conditions. 

Specifically,  the  following  optimal  control  problem  is 
considered:  Minimize  the  functional  I,  which  depends  on  the 

n-vector  state  x(t),  the  m-veetor  control  u(t),  and  the 
p-vector  parameter  n.  The  state  and  the  parameter  are  re- 
quired to  satisfy  r scalar  relations  at  the  initial  point  and 
g scalar  relations  at  the  final  point.  Along  the  interval  of 
integration,  the  state,  the  control,  and  the  parameter  are 
required  to  satisfy  n scalar  differential  equations  and  k 
scalar  nondifferential  relations. 

The  approach  taken  is  a sequence  of  two-phase  cycles, 
composed  of  a conjugate  gradient  phase  and  a restoration 
phase.  The  conjugate  gradient  phase  involves  one  iteration 
and  is  designed  to  decrease  the  value  of  the  functional,  while 
the  constraints  are  satisfied  to  first  order.  During  this  iter- 
ation, the  first  variat ion  of  the  functional  is  minimized,  sub- 
ject to  the  linearized  constraints.  The  minimization  is  per- 
formed over  the  class  of  variations  of  the  control,  the  para- 
meter, and  the  missing  components  of  the  initial  state  which 
are  equidistant  from  some  constant  multiple  of  the  corres- 
ponding variations  of  the  previous  conjugate  gradient  phase. 
The  sequence  of  conjugate  gradient  phases  generated  by  the 
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algorithm  is  such  that,  for  the  special  case  of  a quadratic 
functional  subject  to  linear  constraints,  various  orthogona- 
lity and  conjugacy  conditions  hold.  The  restoration  phase 
involves  one  or  more  iterations  and  is  designed  to  force  cons- 
traint satisfaction  to  a predetermined  accuracy,  while  the 
norm  squared  of  the  variations  of  the  control,  the  parameter, 
and  the  missing  components  of  the  initial  state  is  minimized. 

The  principal  property  of  the  algorithm  is  that  it  pro- 
duces a sequence  of  feasible  suboptimal  solutions:  the  func- 
tions obtained  at  the  end  of  each  cycle  satisfy  the  cons- 
traints to  a predetermined  accuracy.  Therefore,  the  values 
of  the  functional  I corresponding  to  any  two  elements  of  the 
sequence  are  comparable. 

The  stepsize  of  the  conjugate  gradient  phase  is  deter- 
mined by  a one-dimensional  search  on  the  augmented  functional 
J,  while  the  stepsize  of  the  restoration  phase  is  obtained  by 
a one-dimensional  search  on  the  constraint  error  P.  The  con- 
jugate gradient  stepsize  and  the  restoration  stepsize  are 
chosen  so  that  the  restoration  phase  preserves  the  descent 
property  of  the  conjugate  gradient  phase.  Therefore,  the 
value  of  the  functional  I at  the  end  of  any  complete  conjugate 
gradient-restoration  cycle  is  smaller  than  the  value  of  the 
same  functional  at  the  beginning  of  that  cycle.  Of  course, 
restarting  the  algorithm  might  be  occasionally  necessary. 
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is  t=  1.  The  actual  final  time  t,  if  it  is  free,  is  regarded 
as  a component  of  the  vector  parameter  it  to  be  optimized. 

In  this  way,  an  optimal  control  problem  with  variable  final 
time  is  converted  into  an  optimal  control  problem  with  fixed 
final  time. 
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2 . Statement  of  the  Problem 

2.1.  Notation.  Let  t denote  the  independent  variable, 
and  let  x(t),  u(t),  it  denote  the  dependent  variables.  The 
time  t is  a scalar,  the  state  x(t)  is  an  n-vector,  the  con- 

4 

trol  u(t)  is  an  m-vector,  and  the  parameter  it  is  a p-vector. 

The  state  x(t)  is  partitioned  into  vectors  y(t)  and  z(t), 
defined  as  follows:  y(t)  is  an  a-vector  including  those  com- 
ponents of  the  state  that  are  prescribed  at  the  initial  point, 
and  z(t)  is  a b-vector  including  those  components  of  the 
state  that  are  not  prescribed  at  the  initial  point.  Clearly, 
a + b = n. 

2.2.  Optimization  Problem.  With  the  above  notations, 
the  optimization  problem  can  be  stated  as  follows.  Minimize 
the  functional 


= \ f(x,U,TT, 

Jo 


t)  dt  + [h(z,Tr)]Q+  [g(x,n)]1  , 


with  respect  to  the  state  x(t),  the  control  u(t),  and  the 
parameter  tt  which  satisfy  the  differential  constraints 


x - <f>(X,U,TT,t)  = 0, 


0 < t < 1, 


All  vectors  are  column  vectors. 


6 


AAR-145 


the  nondifferential  constraints 

S(x,u,n , t)  = 0,  0 < t < 1,  (3) 

and  the  boundary  conditions 


y (0)  = given, 

(4) 

[u>(z,tt)1=0  , 

(5) 

li|*  (x,H)  ] x = 0 . 

(6) 

In  the  above  equations,  the  quantities  I,f,h,g  are  scalar, 
the  function  tj>  is  an  n-vector,  the  function  S is  a k-vector, 
the  function  at  is  a c-vector,  and  the  function  t|>  is  a q- vector. 
The  number  of  initial  conditions  r = a + c and  the  number  of 
final  conditions  q must  satisfy  the  following  relation: 

r + q < nQ  + n^  + p*  < 2n  + p . (7) 

In  Ineq.  (7) , the  symbol  p*  < p denotes  the  number  of  compo- 
nents of  the  parameter  n present  in  the  boundary  conditions; 
the  symbol  n^  < n denotes  the  number  of  state  variables  pres- 
ent in  the  initial  conditions;  and  the  symbol  n^  < n denotes 
the  number  of  state  variables  present  in  the  final  conditions. 

2.3.  First-Order  Conditions.  From  calculus  of  varia- 
tions, it  can  be  seen  that  the  previous  problem  is  one  of  the 
Bolza  type,  and  it  can  be  recast  as  that  of  minimizing  the 
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5 

agumented  functional 

J=(  [f  + AT(x-4>)+  pTS]dt  + (h  + oTco)  Q + (g  + pTi^)1 

Jo 

= ( (f  - AT4>  + PTS  - ATx)dt+  (-ATx  + h + oTio)  +(\Tx  + g + pTi|/)  , (8) 

Jo 

subject  to  (2)- (6).  Here,  the  n-vector  A(t)  is  a variable 
Lagrange  multiplier,  the  k-vector  p(t)  is  a variable  Lagrange 
multiplier,  the  c-vector  a is  a constant  Lagrange  multiplier, 
and  the  q-vector  p is  a constant  Lagrange  multiplier. 

The  functions  x(t),  u(t),  tt  and  the  multipliers  A(t), 
p(t),  o,  p solving  the  previous  problem  must  satisfy  the 
feasibility  equations  (2) -(6)  and  the  following  optimality 
conditions : 

A - f + <}>  A - S p = 0 , 0<  t<  1,  (9) 

XX  X 


f - 4>  A + s p = o, 

u u uf 


4>  A + S p ) d t + ( h +mo 

it  IT  IT  IT 


0 < t < l , 


JTTVl)l  0' 


( — C + h + io  o ) „ — 0 , 

z z 0 

(*  + gx  + 4'xv« ) 2 = o • 


(10) 

(11) 

(12) 

(13) 


3 The  second  form  of  Eq . (8)  arises  after  the  customary  in- 

tegration by  parts  is  performed. 
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2.4.  Remark . Just  as  the  state  vector  x(t)  is  parti- 


tioned into  an  a-vector  y(t)  and  a b-vector  z(t),  the  multi- 


plier vector  \(t)  is  partitioned  into  an  a-vector  n(t)  and  a 


b-vector  r.(t),  having  the  following  meaning:  n(t)  is  associa- 


ted with  y(t),  and  £(t)  is  associated  with  z(t).  With 


reference  to  Eq.  (12),  C(0)  denotes  the  portion  of  the  initial 


Lagrange  multiplier  A ( 0 ) which  is  associated  with  z(0),  the 


portion  of  the  initial  state  vector  x(0)  which  is  not  pres- 


cribed. 


2.5.  Approximate  Methods.  Since  in  general  the  differ- 


ential system  ( 2 ) — ( 6 ) and.(9)-(13)  is  nonlinear,  approximate 


methods  are  employed  to  find  a solution  iteratively.  In  this 


connection,  let  the  norm  squared  of  a vector  v be  defined  by 


N ( v)  = v v . 


Then,  the  constraint  error  P can  be  written  as 


SN(x-<C)dt+f  N (s  )dt  + N (to)  0 + N (iJj)  x , 

0 Jo 


and  the  error  in  the  optimality  conditions  Q is  given  by 


In  Eq.  (15),  it  is  assumed  that  the  initial  condition  (4)  is 
satisfied . 
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Q = I N(A  " fx  + - Sx 


f1  . fi 

= \ N ( A - f + <b  A - S p ) dt  + 1 N ( f - <j>  A + S p ) dt 
1 X X X 1 u u U 

•'O  •'O 

'[j  (fir  " V + V)dt  + (hn  + V^O  + (gH  + l] 


+ N (-r,  + hz  + wzo)  0 + N (A  + + <J'xVi)  ^ • 


For  the  exact  optimal  solution,  one  must  have 


P = 0V  Q = 0 . 


For  an  approximation  to  the  optimal  solution,  the  following 
relations  are  to  be  satisfied: 


p < , Q < e2  • 


where  and  are  small,  preselected  numbers, 
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3 . Description  of  the  Algorithm 

The  technique  employed  is  characterized  by  a sequence  of 
two-phase  cycles,  composed  of  a conjugate  gradient  phase  and 
a restoration  phase.  These  phases  are  described  below. 

The  conjugate  gradient  phase  is  started  only  when  the 
constraint  error  P satisfies  Ineq.  (18-1).  It  involves  a 
single  iteration,  which  is  designed  to  decrease  the  value  of 
the  functional  I or  the  augmented  functional  J,  while  the 
constraints  are  satisfied  to  first  order.  During  this  itera- 
tion, the  first  variation  of  the  functional  I is  minimized, 
subject  to  the  linearized  constraints.  The  minimization  is 
performed  over  the  class  of  variations  of  the  control  u(t)  and 
the  parameters  n and  z(0)  which  are  equidistant  from  some  con- 
stant multiple  of  the  corresponding  variations  of  the  previous 
conjugate  gradient  phase. ^ 

The  restoration  phase  is  started  only  when  the  constraint 
error  P violates  Ineq.  (18-1) . The  restoration  phase  involves 
one  or  more  iterations.  In  each  restorative  iteration,  the 
objective  is  to  reduce  the  constraint  error  P,  while  the 
constraints  are  satisfied  to  first  order  and  the  norm  squared 


II 


The  sequence  of  conjugate  gradient  phases  generated  by  the 
algorithm  is  such  that,  for  the  special  case  of  a quadratic 
functional  subject  to  linear  constraints,  various  orthogonal- 
ity and  conjugacy  conditions  hold  (see  Section  6 ) . 


■ I ■ ■ .. . <I.i  . ii 
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of  the  variations  of  the  control  u(t)  and  the  parameters  n 
and  z(0)  is  minimized.  The  restoration  phase  is  terminated 
whenever  Ineq.  (18-1)  is  satisfied. 

A complete  conjugate  gradient-restoration  cycle  is  de- 
signed so  that  the  value  of  the  functional  I decreases, while 
the  constraints  are  satisfied  to  the  accuracy  (18-1)  both  at 
the  beginning  and  at  the  end  of  the  cycle.  Finally,  the  al- 
gorithm as  a whole  is  terminated  whenever  Ineqs.  (18)  are 
satisfied  simultaneously. 

3.1.  Remark.  During  each  conjugate  gradient  iteration 
or  restorative  iteration,  use  of  nonlinear  equations  must  be 
avoided.  Therefore,  the  exact  feasibility  equations  ( 2 ) — ( 6 ) 
are  replaced  by  their  corresponding  linearized  approximations. 
These  linearized  approximations  do  not  include  forcing  terms 
in  the  conjugate  gradient  phase,  while  they  do  include  forcing 
terms  in  the  restoration  phase. 

3.2.  Notation.  For  any  iteration  of  the  conjugate 
gradient  phase  or  the  restoration  phase,  the  following  termi- 
nology is  adopted:  x(t),  u(t),  it  denote  the  nominal  functions; 
x(t),  u(t),  it  denote  the  varied  functions;  and  Ax(t),  Au(t), 

Air  denote  the  displacements  leading  from  the  nominal  functions 
to  the  varied  functions.  These  quantities  satisfy  the  relations 


x (t)  = x (t)  + Ax  (t). 


u (t)  = u (t)  + Au  (t). 


fr  = it  + Ait . (19) 


L 
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Because  the  vector  x is  partitioned  into  y and  z,  Eq.  (19-1) 
implies  that 


y 


(t)  = y (t)  + Ay  (t)  , 


z ( t ) =z(t)  + A z ( t ) . 


(20) 


Let  a be  a positive  number  representing  the  stepsize 
\ (either  the  conjugate  gradient  stepsize  or  the  restoration 

stepsize) . Then,  we  define  the  displacements  per  unit  step- 
size  as  follows: 


A ( t)  = Ax  ( t)  /a , B ( t)  = Au  ( t)  /a,  C = Air/ot . ( 21 ) 

The  vector  A is  partitioned  into  vectors  D and  E associated 
with  y and  z,  respectively.  Therefore,  Eq.  (21-1)  implies 
that 

D(t)  = Ay  (t)/«,  E(t)  = Az(t)/o.  (22) 


Upon  combining  (19) -(20)  with  (21) -(22),  we  see  that 

x(t)=x(t)+«A(t),  u(t)=u(t)+aB(t),  Tt  = it  + aC,  (23) 

y(t)=y(t)+oD(t),  z (t)  = z (t)  + «E  (t)  . (24) 


3.3.  Desired  Properties.  The  functions  Ax(t),  Au(t), 

An  must  be  determined  so  as  to  produce  some  desirable  effect 
at  every  iteration,  namely,  the  decrease  of  the  functionals 
I,  and/or  J,  and/or  P.  Thus,  the  following  descent  properties 
are  required: 


u 
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I < I,  and/or  J < J,  and/or  P < P, 


where  1,  J,  P are  associated  with  the  nominal  functions  and 
I,  J,  P are  associated  with  the  varied  functions.  In  turn, 
the  functions  A(t),  B(t),  C are  chosen  so  that 


51 < 0,  and/or  5J  < 0,  and/or  5P  < 0, 


where  the  symbol  5(...)  denotes  the  first  variation.  Then, 
by  choosing  the  stepsize  a sufficiently  small,  the  satisfac- 
tion of  relations  (25)  is  guaranteed.  Ineqs.  (25-1),  (25-2) 

and  (26-1) , (26-2)  characterize  the  conjugate  gradient  phase, 
while  Ineqs.  (25-3)  and  (26-3)  characterize  the  restoration 
phase. 

3.4.  First  Variations.  Next,  we  give  the  expressions 
for  the  first  variations  of  the  functionals  I,  J,  P;  after 
simple  manipulations,  omitted  for  the  sake  of  brevity,  they 
take  the  form®'** 


SI 

T t T T T T T 

(f  A + f B + f C)dt  + (h  E + h C) + (q  A + g C)  . , (27) 

x u it  z it  0 ^x  it  1 

0 


implicit  in  Eqs.  (27)  — (29)  is  the  assumption  D ( 0 ) =0. 

*The  first  variation  of  the  augmented  functional  J is  computed 
by  varying  the  functions  x(t),  u(t),  it,  while  holding  the 
multipliers  X(t),  o(t),  o,  |i  unchanged. 
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and 


,SJ/a=  l (-X  + f - <fc\  + S vp)TAdt  + J1  (f  - <p  \ + S p)TBdt 

JQ  XX  X JqUU  U 

+ (10  <f™-*»X  + sn*'>dtt  <h„+V’>0*  (9„  +*„">1]TC 

* [(-C  + h2  +<^zo)Te]o+  [u  + 9x  + V)Ta]1  ' 


: 28 ) 


and 


• T . T T T f 1 rp  rr>  m m 

(SP/2ct  = l (x-t)1  (A~4'xA-4'uB-l)'ijC)dt  + l S (S„A+S,  B+S„C)dt 

J0 


= j1  <x-< 
'0 


X U T1 


+ £u>T(w^E  + W^C)  ] 


[tT(^A+g^C)]i 


(29) 


3.5.  Remark.  For  the  purpose  of  this  report,  Eqs. 

(27)- (29)  must  be  completed  by  one  of  the  following  relations 


■2=  (‘ 


K/a"  = |oBTBdt  + CTC+  (ETE)0, 


(30) 


or 


K/oi2  = ( B->  B)  T ( B- y B)  dt+  (C- yC)  T (C- yC)  + j^(E- yE)  T ( E- yE) J t (31) 


i 


u 

[J 

j 


.) 


I J 


J 


- 


u • 
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which  constitute  a measure  of  the  overall  change  of  the  con- 
trol u(t)  and  the  parameters  it  and  z(0).  Equation  (30)  is 
employed  in  the  restoration  phase,  and  Eq.  (31)  is  employed 
in  the  conjugate  gradient  phase.  In  Eq.  (31) , the  functions 
A ( t ) , B ( t ) , C pertain  to  the  present  conjugate  gradient  phase, 

AAA 

and  the  functions  A(t),  B(t),  C pertain  to  the  previous  con- 
jugate gradient  phase.  The  symbol  y denotes  a scalar,  non- 
negative quantity,  called  the  directional  coefficient.  Its 
specification  is  discussed  in  Section  6 . 
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4 . Restoration  Phase 

4.1.  Linearized  Equations.  Let  x(t),  u(t),  it  denote 
nominal  functions  satisfying  (4),  but  not  necessarily  (2) -(3) 
and  (5) -(6).  To  first  order,  the  perturbations  per  unit 
stepsize  A(t) , B(t),  C must  satisfy  the  linearized  constraint 
equations 

A - 0^A  - <|>^B  - (p^C  + (x  - <f>)  = 0,  °<  t<  1,  (32) 


i 


E^  + SJB  + S^C  + S = 0,  0<t<l, 


(33) 


D(0)  = 0, 


(34) 


( U)Je  + U)Jc  + W ) g = 0 , 


(35) 


(t^A  + tJ^C  + tM^  0.  (36) 

4.2.  Descent  Property.  The  linearized  equations 
(32)- (36)  admit  an  infinite  number  of  solutions,  each  of 
which  is  characterized  by  a descent  property  in  the  constraint 
error  P.  This  descent  property  can  be  seen  by  combining  (29) 
with  (32)— (36) c the  first  variation  of  P becomes 

<$P  = -2otP . (37) 


L 


L I 


) 


■ 


J 
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Since  P>0,  Eq.  (37)  shows  that  6P<0;  hence,  for  a suffi- 
ciently small,  a decrease  in  the  constraint  error  P is  guar- 


anteed. 


4.3.  Special  Variations.  Among  the  infinite  number  of 
solutions  of  Eqs.  (32)- (36),  the  one  that  minimizes  the 
functional  (30)  is  selected.  Thus,  we  seek  the  minimum  of 
the  quadratic  functional  (30),  with  respect  to  the  functions 
A ( t ) , B ( t ) , C which  satisfy  the  linearized  constraints 
(32) - (36) . 

By  applying  standard  techniques  of  optimal  control  theory 
or  calculus  of  variations,*  the  following  optimality  conditions 


are  obtained: 


B=  V ‘ V' 


0 < t < 1 , 


C=  L(V"V)dt_  (a)7Ta)0‘  (VJ)1' 


E(0)  = (£  - U)zO)  g , 


A + (|>XA  - S^p  = 0 , 0 < t < 1, 


(X  + i|>xu)  i = 0 


Summarizing,  we  seek  functions  A(t),  B(t),  C and  multipliers 
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X (t)  , p(t),  o,  (i  which  satisfy  the  linearized  constraints 
(32) -(36)  and  the  optimality  conditions  (38) -(42). 

4.4.  Linear,  Two-Point  Boundary-Value  Problem.  The 
technique  used  to  solve  the  LTP-BVP  (32)  — (36)  and  (38)  — (42) 
is  a forward  integration  scheme  in  combination  with  the  method 
of  particular  solutions  (Refs.  6-8).  The  technique  requires 
the  execution  of  n + p + 1 independent  sweeps  of  the  differen- 
tial system  (32)  — (36)  and  (38)  — (42) , each  characterized  by  a 
different  value  of  the  (n  + p)-vector  w,  whose  components  are 
the  n components  of  the  initial  multiplier  \(0)  and  the  p 
components  of  the  parameter  C. 

The  generic  sweep  is  started  by  assigning  particular 
values  to  the  components  of  w,  that  is,  the  components  of  the 
vectors  X (0)  and  C.  With  X (0)  known,  C ( 0 ) is  known.  There- 
fore, Eqs.  (35)  and  (40)  constitute  a system  of  b + c linear 
relations  in  which  the  unknowns  are  the  b + c components  of 
the  vectors  E(0)  and  o.  For  this  system  to  have  a unique 
solution,  the  following  disequation  must  hold:^ 


^Disequation  (43)  is  obtained  from  (35)  and  (40)  after  elimin- 
ation of  E (0) . The  resulting  linear  equation  in  a admits  a 
unique  solution  providing  (43)  is  satisfied. 


r-*i 
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With  E ( 0 ) known  and  because  of  (34),  the  vector  A(0)  is  known. 
Then,  A(t)  and  X (t)  together  with  B(t)  and  p(t)  are  obtained 
by  forward  integration  of  (32)  and  (41)  , subject  to  (33)  and 
(38).  Note  that,  at  each  time  station  t,  0<t<l,  Eqs.  (33) 
and  (38)  constitute  a system  of  m + k linear  relations  in  which 
the  unknowns  are  the  m + k components  of  the  vectors  B(t)  and 
p(t).  For  this  system  to  have  a unique  solution,  the  follow- 
ing disequation  must  hold:"*"^ 


det 


[susu]''°- 


0 < t < 1 , 


(44) 


i. 

L 

V 

[ 

[ 


As  a result  of  the  procedure,  the  sweep  is  completed:  for  the 
arbitrary  value  assigned  to  w,  it  leads  to  the  satisfaction  of 
all  of  the  equations  of  the  system  (32) -(36)  and  (38) -(42), 
except  Eqs.  (36),  (39),  (42). 

In  order  to  satisfy  Eqs.  (36)  , (39) , (42)  and  because 

the  system  (32) -(36)  and  (38) -(42)  is  nonhomogeneous , n + p+1 
independent  sweeps  must  be  executed  employing  n + p+1  differ- 
ent vectors  w^,  i=l,...,n  + p+l.  The  first  n + p sweeps  are 

performed  by  choosing  the  vectors  w^, ,wn+p  to  *"*ie 

columns  of  the  identity  matrix  of  order  n + p.  The  last  sweep 


11 


Disequation  (44)  is  obtained  from  (33)  and  (38)  after  elim- 
ination of  B ( t) . The  resulting  linear  equation  in  p(t) 
admits  a unique  solution  providing  (44)  is  satisfied. 
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is  executed  by  choosing  wn  + p + to  be  the  null  vector.  As 
a result,  one  generates  the  functions  and  multipliers 

Ai  (t),  Bi  (t),  Ci,Xi  (t),  p±  (t),  o i , i = 1, n + p + 1 . (45) 

Now,  we  introduce  the  n + p+1  undetermined,  scalar  con- 
stants k^  and  form  the  linear  combinations 

A(t)  = T.kiAi  (t)  , B(t)  = EkiBi  (t)  , C = Ek^  , (46) 

X(t)  = Zki\i(t)  , p(t)  = Hkipi(t)  , o = Ekiai,  (47) 

where  the  summations  are  taken  over  the  index  i.  The  n + p+1 
coefficients  k^  and  the  q components  of  the  multiplier  \i  are 
obtained  by  forcing  the  linear  combinations  (46) -(47)  to  sat- 
isfy Eqs.  (36) , (39) , (42) , together  with  the  normalization 
condition  (Ref.  6) 


Iki  = 1 . (48) 

Once  the  constants  k^  are  known,  the  solution  of  the  LTP-BVP 
(32)— (36)  and  (38)-(42)  is  given  by  (46)-(47). 

4.5.  Restoration  Stepsize.  With  the  functions  A(t), 

B ( t) , C known,  the  one-parameter  family  of  varied  functions 
(23)  can  be  formed.  For  this  one-parameter  family,  the  con- 
straint error  (15)  becomes  a function  of  the  form 


iJ 


) 

0 


L 

1. 

L 

li 

0 

li 
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Then,  the  stepsize  a must  be  selected  so  that  the  following 
relations  are  satisfied: 

P (a)  < P{0) , i (a)  > 0.  (50) 

Satisfaction  of  Ineq.  (50-1)  is  possible  because  of  the  des- 
cent property  of  the  restoration  phase.  Ineq.  (50-2)  is 
required  for  problems  with  free  final  time. 

In  order  to  achieve  satisfaction  of  (50),  a bisection 
process  is  applied  to  the  restoration  stepsize  a,  starting 
from  the  reference  stepsize  u =1.  This  reference  stepsize 
has  the  property  of  yielding  one-step  restoration  for  the  case 
where  the  constraints  (2) -(6)  are  linear. 

4.6.  Iterative  Procedure  for  the  Restoration  Phase. 

The  descent  property  (37)  of  the  restoration  phase  guarantees 
satisfaction  of  Ineq.  (50-1)  at  the  end  of  any  iteration,  but 
not  satisfaction  of  Ineq.  (18-1).  Therefore,  the  restoration 
algorithm  must  be  employed  iteratively  until  Tneq.  (18-1)  is 
satisfied.  At  this  point,  the  restoration  phase  is  termina- 
ted. 


i 


, 
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5 . Conjugate  Gradient  Phase:  General  Case 

5.1.  Linearized  Equations.  Suppose  that  nominal  func- 
tions x(t),  u(t),  7i  are  available,  which  satisfy  (2)  — (6)  . To 
first  order,  the  perturbations  per  unit  stepsize  A(t),  B(t), 

C must  satisfy  the  linearized  constraint  equations 


A- <j>^A-  <|>Jb-  <i>^c=  °,  0<  t<  1,  (51) 

STA  + STB  + STC  = 0,  0<t<l,  (52) 

X U 77  ' - 

D ( 0 ) = 0,  (53) 

(u^E  + u)^C)0=  0,  (54) 

(^A+^C)1  = 0.  (55) 


5.2.  Special  Variations.  Among  the  infinite  number  of 
solutions  of  Eqs . (51)-(55),  the  one  that  minimizes  the  func- 

tional (27)  is  selected.  Thus,  we  seek  the  minimum  of  the 
linear  functional  (27) , with  respect  to  the  functions  A(t) , 

B ( t ) , C which  satisfy  the  linearized  constraints  (51)— (55) 
and  the  quadratic  isoperimetric  constraint  (31). 

By  applying  standard  techniques  of  optimal  control  theory 
or  calculus  of  variations,  the  following  optimality  conditions 


are  obtained: 


12 


B = yB  - (f  — <j>  A + S p ) , 0<t<l, 

u u u ' - 


(56) 


C = yC 


-[£«.- 


♦ ^A  + S^p ) dt  + (h^  + 0)^0)  Q + (g 


„ + V>i]' 


(57) 


E(0)  = yE(0)—  (-  c-*-  hz+  a)za ) 0 , 


(58) 


A - f x + 4>xA  - Sxp  = 0 , 0 < t < 1 , 


(59) 


(A  + gx  + i|)xm  ) x = 0 . 


(60) 


Summarizing,  for  a given  value  of  the  directional  coefficient 
Y,  we  seek  functions  A(t) , B(t),  C and  multipliers  A(t),  p(t), 
a,  y which  satisfy  the  linearized  constraints  (51)- (55)  and 
the  optimality  conditions  (56) -(60). 

5.3.  Isoperimetric  Constraint.  In  the  light  of  (56)-(60), 
the  error  in  the  optimality  conditions  (16)  reduces  to 

Q=  j (B-YB)T(B-YB)dt  + (C-yC)T(C-yC)  + [(E-yE)T(E-yE)]q  . (61) 


12 


In  Eqs . (56)— (60) , the  multiplier  v associated  with  the 
isoperimetric  constraint  (31)  is  set  at  the  level  2v=l. 


* 


Consequently,  the  following  relation  ties  the  isoperimetric 
constant,  the  stepsize,  and  the  error  in  the  optimality  con- 
ditions : 

K = a2Q  . (62) 

Clearly,  to  assign  values  to  the  isoperimetric  constant 
is  the  same  as  assigning  values  to  the  stepsize.  However, 
there  is  no  clear-cut  way  of  determining  a priori  convenient 
values  for  the  isoperimetric  constant.  Therefore,  the  imple- 
mentation of  the  conjugate  gradient  algorithm  becomes  simpler 
if  one  avoids  evaluating  \t  in  terms  of  K through  (62)  and 
assigns  values  to  a directly. 


5.4.  Descent  Property.  When  the  variations  defined  by 
(51) -(60)  are  employed,  the  first  variation  of  the  augmented 
functional  (28)  becomes 

5j  = -a (Q  + yZ) , (63) 

where  Q is  given  by  (61)  and 

Z=\  (B  - yB)TBdt  + (C  - yC)TC  + [(E  - yE)TE  ]0  • (64) 

>0 

For  the  first  iteration  of  the  conjugate  gradient  phase, 
one  sets 


Y=0,  (65) 


I 
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with  the  implication  that 


6J  - -aQ  . ( 66 ) 

Since  Q N 0 , Eq.  (66)  shows  that  OvT  < 0.  Hence,  for  a suffi- 
ciently small,  it  is  guaranteed  that  the  augmented  functional 
J decreases. 

For  subsequent  iterations,  one  sets  ^ *0.  More  specifi- 
cally, the  directional  coefficient  must  be  such  that 

i > 0,  (67) 

and  its  proper  value  is  discussed  in  Section  6.  At  any  rate, 
Eq.  (63)  shows  that  <$  j v 0 providing 

0 + i 7,  > 0 . ((,«) 

Hence,  for  a sufficiently  small,  it  is  guaranteed  that  the 
augmented  functional  J decreases  as  long  as  Itieq,  (68)  is 
satisfied.  If  Ineq.  (68)  is  violated,  the  descent  property 
on  J no  longer  holds,  and  the  conjugate  gradient  phase  must 
be  restarted  by  resetting  the  directional  coefficient  ^ at  the 
level  (65) . 

5.5.  Linear,  Two-Point  Boundary-Value  Problem.  For  a 
given  value  of  the  directional  coefficient  ^ , the  technique 
used  to  solve  the  LTP-BVP  (51) -(60),  associated  with  the  con- 
jugate gradient  phase,  is  analogous  to  that  described  for  the 


I 


26 


AAR-145 


restoration  phase  (see  Section  4.4);  hence,  it  is  not  re- 
peated, for  the  sake  of  brevity. 

5.6.  General  Solution.  Next,  assume  that  two  particu- 
lar values  are  given  to  the  directional  coefficient  y,  for 
instance , 


Y * = 0 and  Y **  = 1 . (69) 

Denote  by 

A*(t),  B*(t),  C*,  X* (t) , P*(t),  a*  ,u*  (70) 

and 

A**  (t),  b**(t),  C**,X**(t),  P**(t),  o**,u**  (71) 

the  particular  solutions  of  the  LTF-BVP  (51) -(60)  correspond- 
ing to  (69-1)  and  (69-2)  , respectively.  Simple  manipulations, 
omitted  for  the  sake  of  brevity,  show  that  the  general  solu- 
tion of  (51) -(60),  valid  for  any  value  of  the  directional 
coefficient  y,  can  be  written  as 

A ( t)  = A*(t)  + Y(A**(t)  - A* ( t ) ) , (72-1) 


B ( t ) = B* (t)  + Y lB**(t)  - B* (t) ] , (72-2) 


C = + Y (C**  - CA)  , 


(72-3) 
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A(t)  = A*(t)  + YU**(t)  - A*  ( t)  ] , 


(73-1) 


p(t)  = p* (t)  + Y tP** (t)  - p* (t) ] , 


(73-2) 


° = °*  + 7(o**  - a*)  , 


(73-3) 


u = u*  + Y(y**-y*) 


(73-4) 


As  a conclusion,  the  general  solution  of  (51) -(60)  requires 
that  two  sets  of  n + p + 1 sweeps  be  executed,  one  leading  to 
the  particular  solution  (70)  and  one  leading  to  the  particu- 
lar solution  (71) . 

5.7.  Stepsize  and  Directional  Coefficient.  With  the 
functions  A(t),  B(t),  C known,  the  following  two-parameter 
family  of  varied  functions  can  be  formed: 

x(t)  = x(t)  +a{A*(t)  +Y[A**(t)  - A*(t))j-  , (74-1) 

u(t)  =u(t)  +ot|B*(t)  +YtB**(t)  — B*(t)]}  , (74-2) 


if  = 7T  + a [C*  + y (C**  - C*)  ] . 


(74-3) 


On  the  other  hand,  the  multipliers  A ( t) , p (t) , a , p form  the 
one-parameter  family  (73) . Upon  using  (73)  and  (74)  , we  see 
that  the  augmented  functional  (8)  takes  the  form 
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J « J (a , y)  . (75) 

Therefore,  the  optimum  values  of  a and  y satisfy  the  relations 
Ja(a,Y)  = 0,  (a , y ) = 0 . (76) 

Since  the  simultaneous  determination  of  a and  > might  be 
be  expensive  computationally,  we  proceed  in  a different  way. 
First,  we  determine  an  approximate  value  of  the  directional 
coefficient  > , based  on  the  consideration  of  the  linear- 
quadratic  model  (Section  6).  Once  y is  known,  the  two- 
parameter  family  (75)  reduces  to  the  one-parameter  family 

J = J(a)  . (77) 

Then,  the  optimum  stepsize  a satisfies  the  relation 

J (a)  = 0 , (78) 

whose  numerical  solution  can  be  obtained  using  quadratic  in- 
terpolation or  cubic  interpolation  (Ref.  9) . 
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6 . Conjugate  Gradient  Phase:  Linear-Quadratic  Case 

In  the  previous  section,  we  analyzed  the  conjugate 
gradient  phase  in  general,  regardless  of  the  analytical  form 
of  the  functional  (1)  and  the  constraints  (2)-(6).  In  this 
section,  we  consider  the  linear-quadratic  case,  that  is,  the 
case  where  the  functional  (1)  is  quadratic  and  the  constraints 
(2)- (6)  are  linear. 

6.1.  General  Solution.  Under  the  assumption  of  linear 

constraints,  it  can  be  verified  that  the  particular  solutions 

(70)  and  (71)  satisfy  the  relations 

« 

A**  (t)-A*  (t)=  A(t),  B**(t)-B*(t)=B(t),  C**-C*=C,  (79) 

A**  (t)-A*  (t)  = 0»  P**(t)~P*(t)  = 0,  cr**-a*  = 0,  P**“P*=0.  (80) 

As  a consequence,  Eqs.  (72) -(73)  reduce  to 


A(t)  = A*  (t)  + YA(t)  , B(t)=  B*  (t)  + yB(t),  C = C*+  yC  , (81) 

A(t)=A*(t),  p(t)=p*(t),  a = °*  / P = M*  . (82) 

This  means  that  the  general  solution  of  the  LTP-BVP  (51) -(60) 
can  be  obtained  by  executing  only  one  set  of  n + p + 1 sweeps, 
namely,  the  set  of  sweeps  leading  to  the  solution  (70) . By 
the  way,  this  is  the  solution  corresponding  to  (69-1)  , namely, 
the  solution  associated  with  the  ordinary  gradient  phase  of  Ref . 5. 
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6.2.  Local  Orthogonality  Conditions.  Under  the  assump- 
tion of  linear  constraints,  the  two-parameter  family  (74) 
simplifies  to 


x(t)  = x(t)  + a[A*(t)  + yA(t)  ] , 


(83-1) 


u(t)  — u ( t)  +a[B*(t)  + yB(t)], 


(83-2) 


ir  = it  + a (C*  + YC)  . 


(83-3) 


On  the  other  hand,  the  multipliers  X(t),  p(t),  a,  u are  given 
by  Eqs . (82).  Upon  using  (82)  and  (83),  we  see  that  the 
augmented  functional  (8)  still  takes  the  form  (75).  Hence, 
the  optimum  values  of  a and  y still  satisfy  the  relations  (76). 

After  laborious  manipulations,  omitted  for  the  sake  of 
brevity,  Eqs.  (76)  lead  to  the  following  local  orthogonality 
conditions : 


J B*Bdt  + C*C  + (E*E)  0=0, 


l B*Bdt  + C 

Jn 


cic+  (e;e)0-  o, 


(84-1) 


(84-2) 


with  the  implication  that 


®*B*dt  + C*C*+  (E*E*)  q=0, 

"ft 


(84-3) 
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Here,  the  adjective  local  is  employed  to  mean  that  Eqs.  (84) 
involve  vectors  B(t),  C,  E(0)  which  are  solutions  of  (51) -(60) 
computed  for  the  present  iteration  and  the  previous  iteration; 
they  also  involve  vectors  B*(t),  C*,  E*(0)  which  are  solutions 
of  (51)  — (60)  for  y = 0 computed  for  the  present  iteration  and 
the  next  iteration. 

6.3.  Local  Conjugacy  Condition.  Let  w(t)  and  M(t) 
denote  the  vectors 


w ( t ) = u ( t ) 


A ( t ) 

M(t)  = B ( t ) 


Let  f , g , h denote  the  Hessian  matrices  of  the  functions 
ww  ww  ww 

f,g,h  with  respect  to  the  vector  w.  With  this  notation,  and 
under  the  assumption  of  linear  constraints  and  quadratic 
functional,  Eqs.  (76)  lead  to  the  following  local  conjugacy 


condition: 


[ MTf  Mdt+(MTh  M)  + (MTg  M)  =0. 
\ ww  ww  0 5ww  1 

Jo 


Here,  the  adjective  local  is  employed  to  mean  that  Eq.  (86) 
involves  vectors  M(t),  that  is,  vectors  A(t),  B(t),  C, which 
are  solutions  of  (51) -(60)  computed  for  the  present  iteration 
and  the  previous  iteration. 


» 
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6.4.  Stepsize  and  Directional  Coefficient.  After 
observing  that 

M ( t)  = M* (t)  + yM(t) , (87) 

and  after  accounting  for  Eqs.  (84)  and  (36),  it  can  be  shown 
that  the  optimal  values  of  y and  a are  given  by 

Y = Q/Q  , a = Q/R  , (88) 


where 

(1 

Q = )0  B*B*dt  + C*C*+(E*E*) 0 , (89-1) 

Q=  B*B*dt  + C*C*  + (E^E*)  Q , (89-2) 

R=  SoMTfwwMdt+  (MThwwM)0  + (MT5wwM)l*  (89"3) 

Clearly,  the  optimal  directional  coefficient  y is  the 
ratio  of  the  error  in  the  optimality  conditions  Q for  the 
present  conjugate  gradient  iteration  to  the  error  in  the 

A 

optimality  conditions  Q for  the  previous  conjugate  gradient 
iteration.  These  quantities  are  known,  since  they  involve 
vectors  B* (t) , C*,  E*(0)  which  are  solutions  of  (51)-(60)  for 
y = 0 computed  for  the  present  iteration  and  the  previous 
iteration 


» 


I. 

r. 
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t . 


L 

L 


C 

0 


With  y known,  the  vector  M(t)  is  computed  with  (87)  and 
the  scalar  quantity  R is  computed  with  (89-3) . Then,  the 
optimal  stepsize  a is  evaluated  with  (88-2) . Clearly,  the 
optimal  stepsize  a is  the  ratio  of  the  error  in  the  optimal- 
ity conditions  Q to  the  scalar  quantity  R,  which  constitutes 
a measure  of  the  curvature  of  the  functional  (1) . Both  Q 
and  R are  computed  employing  quantities  pertaining  to  the 
present  conjugate  gradient  iteration. 

6.5.  Descent  Property.  Under  the  assumption  of  linear 
constraints,  Eq.  (64)  simplifies  to 

T'' 

Z = ^ B;Bdt  + c*c  + (e;e) 0 . (90) 

Because  of  the  local  orthogonality  condition  (84-1)  written 
for  the  previous  iteration,  Eq.  (90)  yields 

Z=0.  (91) 

As  a consequence,  Eq.  (63)  reduces  to 

SJ=-aQ,  (92) 

where  the  error  in  the  optimality  conditions  Q is  given  by  Eq. 

(89-1) . Equation  (92)  holds  for  any  conjugate  gradient  iter- 
ation and  shows  that,  since  Q > 0,  we  have  SJ  < 0.  Hence,  for 
a sufficiently  small,  it  is  guaranteed  that  the  augmented 
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functional  0 decreases.  In  conclusion,  for  the  linear- 
quadratic  case,  the  restart  procedure  mentioned  in  Section 
5.4  never  occurs.  This  means  that  the  directional  coeffi- 
cient 'i  is  set  at  the  level  (65)  only  for  the  first  conjugate 
gradient  iteration. 

6.6.  General  Orthogonal i ty  and  Conjugacy  Conditions . 

Now,  assume  that  the  algorithm  described  by  Kgs.  (51) -(60) 
and  (83)  is  employed,  starting  with  some  feasible  nominal 
functions.  Further,  assume  that  the  directional  coefficient 
is  set  at  the  level  (65)  for  the  first  conjugate  gradient 
iteration  and  at  the  level  (88-1)  for  any  subsequent  conjugate 
gradient  iteration.  Under  thos  ' assumptions  and  for  the 
linear-quadratic  case,  one  can  generalize  the  local  orthogo- 
nality conditions  (84)  and  the  local  conjugacy  condition 
(86)  as  follows: 


m mi  nl 

B n dt  + crc  + (k  k ) n = o, 

* p * p * p 0 


(93-1) 


T T T 

B B dt  + C C f (E  E ) „ 0, 

* *p  * *p  * *p  0 


(93-2) 


and 


[l  MTf  M dt  + (MTh  M )„  + (MTg  M ),  =0, 
ww  p ww  p 0 ww  p 1 


(93-3) 
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where  the  subscript  p denotes  any  conjugate  gradient  itera- 
tion preceding  the  present  conjugate  gradient  iteration. 

While  these  equations  do  not  guarantee  convergence  in  a finite 
number  of  steps,  they  do  guarantee  that  the  algorithm  gener- 
ates a sequence  of  linearly  independent  vectors  M(t),  that 
is,  a sequence  of  linearly  independent  variations  per  unit 
stepsize  A(t),  B(t),  C. 
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7 . Conjugate  Gradient  Phase:  Practical  Implementation 

In  this  section,  we  summarize  the  results  of  Sections 
5-6  and  suggest  practical  ways  of  utilizing  these  results, 
within  the  following  operating  rules:  (i)  the  use  of  second 
derivatives  must  be  avoided;  and  (ii)  the  simultaneous  de- 
termination of  the  directional  coefficient  y and  the  stepsize 
ct  must  be  avoided.  We  refer  to  the  general  case  where  the 
functional  (1)  is  nonquadratic  and/or  the  constraints  (2) -(6) 
are  nonlinear. 

7.1.  Auxiliary  Functions.  The  first  step  is  to  solve 
Eqs . (51) -(60)  for  a fictitious  value  of  the  directional 

coefficient,  namely. 


Y*  = 0 . (94) 

Using  the  solution  technique  of  Section  4.4,  we  obtain  the 
following  auxiliary  functions  and  multipliers: 

A*  (t),  B*(t),  C*,A*(t),  P*(t),  o*,  y*  , (95) 

7.2.  Directional  Coefficient.  The  second  step  is  to 
compute  the  actual  value  of  the  directional  coefficient  y . 
For  the  first  conjugate  gradient  phase,  we  set 

Y = 0.  (96) 

For  subsequent  conjugate  gradient  phases,  we  set 
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Y = Q/Q  t 


(97) 


where 


Q = 

( 1 m ip 

\ BXdt  + C*C*  + 

Jo 

(eJe*)0, 

(98-1) 

| 

[ J aTa  'T' 

Q = 1 

1 BjB*dt  + C*C*  + 

J0 

(e;e*)0  . 

(98-2) 

In  Eqs.  (97)  — (98) r the  symbols  Q and  Q denote  the  errors  in 
the  optimality  conditions  for  the  present  conjugate  gradient 
phase  and  the  previous  conjugate  gradient  phase,  respectively. 
The  directional  coefficient  (97)  is  acceptable  only  if 


Ja(0)  = - (Q  + YZ)  < 0 , 


(99) 


where  Q is  given  by  (98-1)  and  Z is  given  by 


= \1bIb 

J0 


Bdt  + C*C  + (E*E)  . 


(100! 


If  Ineq.  (99)  is  violated,  then  the  directional  coefficient 
(97)  must  be  discarded  and  replaced  by  the  value  (96) . This 
means  that  the  algorithm  must  be  restarted  by  replacing  the 
conjugate  gradient  phase  with  an  ordinary  gradient  phase. 

7.3.  Basic  Functions.  The  third  step  is  to  compute  the 
basic  functions  A(t),  B(t),  C and  the  multipliers  \(t),  p(t), 
a,  u • This  is  done  with  the  following  formulas: 
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A(t)  = A* (t)  + YA(t),  B(t)  = B*  (t)  + ')B(t),  C=C*+YC,  (101) 

A (t)  = A*  (t),  P(t)=p*(t),  o = o*,  u = U*.  (102) 

Therefore,  in  the  practical  implementation  of  the  algorithm, 
the  basic  functions  A (t),  B(t),  C are  computed  using  the  formu- 
las derived  under  the  assumption  of  linear  constraints. 

7.4.  Stepsize . With  the  basic  functions  (101)  known, 
we  consider  the  one-parameter  family  of  varied  functions 

x(t)=x(t)+aA(t),  u(t)=u(t)  + aB(t),  it  = n + aC . (103) 

After  substitution  of  Eqs.  (102)-(103)  into  (8)  and  (15),  the 
following  functions  of  the  stepsize  are  obtained: 

J = J (ct)  , P = P ( a ) . (104) 

Then,  a one  dimensional  search  scheme  is  applied  to 
(104-1),  and  a value  of  the  stepsize  a is  selected  for  which 
the  following  relations  are  satisfied: 

J (a)  < J (0)  , P (a)  < P*  , f ( a ) > 0 , (105) 

where  t is  the  final  time  and  P*  is  a preselected  number, 
not  necessarily  small.  Satisfaction  of  Ineq.  (105-1)  is  pos- 
sible because  of  the  descent  property  of  the  conjugate 
gradient  phase.  Ineq.  (105-2)  is  introduced  to  prevent 
excessive  constraint  violation.  And  Ineq.  (105-3)  is  required 
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[ 

c 

[ 

i 


for  problems  with  free  final  time. 

Prior  to  the  satisfaction  of  (105) , a scanning  process 

is  employed,  leading  to  the  bracketing  of  the  minimum  point 

for  J(a).  This  operation  is  then  followed  by  a Hermitian 

cubic  interpolation  process  (Ref.  9) , which  is  stopped  whenever 

1 3 

the  following  relation  is  satisfied: 


r 

L 


E 


i: 

c 

c 


|Ja(a)|<c3  or  | Ja(«)/Jrt(0) | < e4,  (106) 

subject  to  an  upper  limit  for  the  number  of  search  steps  N . 

s 

Once  a stepsize  « has  been  selected  consistently  with  either 
(106)  or  the  prescribed  upper  limit  for  the  number  of  search 
steps,  Ineqs.  (105)  must  be  checked.  If  satisfaction  occurs, 
then  the  stepsize  ci  is  accepted.  If  any  violation  occurs, 
then  the  stepsize  aQ  must  be  bisected  progressively  until 
satisfaction  of  (105)  is  finally  achieved. 


13 


The  symbols  i ^ and  t ^ denote  small,  preselected  numbers, 


1.1  It.  VtVojI 
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8 . Descent  Property  of  a Cycle 

A descent  property  exists  for  a complete  conjugate 

gradient-restoration  cycle  under  the  assumption  of  small 

stepsizes.  Let  denote  the  conjugate  gradient  stepsize  and 

otr  the  restoration  stepsize.  Simple  manipulations,  omitted 

for  the  sake  of  brevity,  show  that  the  conjugate  gradient 

corrections  are  of  CHcx^),  while  the  restoration  corrections 
2 

are  of  0(a  a ).  Hence,  for  a sufficiently  small,  the  restor- 
r g g 

ation  corrections  are  negligible  with  respect  to  the  conjugate 
gradient  corrections.  Therefore,  the  restoration  phase  pre- 
serves the  descent  property  of  the  conjugate  gradient  phase. 

More  specifically,  let  I ^ , I0,  1^  denote  the  values  of 
the  functional  (1)  at  the  beginning  of  the  conjugate  gradient 
phase,  at  the  end  of  the  conjugate  gradient  phase,  and  at  the 
end  of  the  subsequent  restoration  phase.  Note  that  1^  and  I0 
are  not  comparable,  since  the  constraints  are  not  satisfied 
to  the  same  accuracy.  On  the  other  hand,  I3  and  1^  are  com- 
parable, and  the  conjugate  gradient  stepsize  u can  be  selected 
so  that 


I3  < Ix  . (107) 

This  inequality  constitutes  the  descent  property  of  a complete 
conjugate  gradient-restoration  cycle.  In  order  to  enforce  it, 
one  proceeds  as  follows.  At  the  end  of  the  restoration  phase, 


■ — 


one  must  verify  Ineq.  (107).  If  it  is  satisfied,  the  next 
conjugate  gradient  phase  is  started;  otherwise,  the  previous 
conjugate  gradient  stepsize  is  bisected  as  many  times  as 
needed  until,  after  restoration,  Ineq.  (107)  is  satisfied. 
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9 . Summary  of  the  Algorithm 

The  sequential  conjugate  gradient-restoration  algorithm 
solves  optimal  control  problems  involving  a functional,  sub- 
ject to  differential  constraints,  nondifferential  constraints, 
and  general  boundary  conditions.  The  algorithm  is  composed 
of  a sequence  of  cycles,  each  cycle  consisting  of  two  phases, 
a conjugate  gradient  phase  and  a restoration  phase.  The 
objective  of  each  cycle  is  to  decrease  the  functional  I, 
while  the  constraints  are  satisfied  to  the  predetermined 
accuracy  (18-1). 

The  decision  parameters  controlling  the  algorithm  are 
the  constraint  error  P and  the  optimality  condition  error  Q 
(see  Eqs . (15)  and  (16) J.  If  P violates  Ineq.  (18-1),  the 

algorithm  executes  a restoration  phase.  If  P satisfies 
Ineq.  (18-1)  and  Q violates  Ineq.  (18-2),  the  algorithm  executes  a 
conjugate  gradient  phase.  Finally,  if  P and  Q satisfy  Ineqs. 
(18),  the  algorithm  stops:  convergence  has  been  achieved. 

9.1.  Restoration  Phase.  This  phase  involves  one  or 
more  iterations  and  can  be  summarized  as  follows. 

(a)  Assume  nominal  functions  x(t),  u(t),  tt  which  satis- 
fy condition  (4),  but  violate  at  least  one  of  conditions 
(2) - (3)  and  (5) - (6)  . 


' 


1 


J 


(b)  For  the  nominal  functions,  solve  the  LTP-BVP 
(32)- (36)  and  (38)- (42)  using  the  method  of  particular 
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solutions.  In  this  way,  obtain  the  functions  A(t),  B(t), 
and  the  multipliers  A(t),  p(t),  o , p • 

(c)  Using  the  functions  in  (b) , compute  the  restoration 
stepsize  by  a one-dimensional  search  on  the  constraint  error 
P (a)  . To  this  effect,  perform  a bisection  process  on  a, 
starting  from  a = 1 , until  Ineqs.  (50)  are  satisfied. 

(d)  Once  the  restoration  stepsize  is  known,  compute  the 

varied  functions  x(t),  u(t),  it  with  Eqs . (23). 

(e)  Verify  whether  the  varied  functions  in  (d)  satisfy 
Ineq.  (18-1).  If  this  is  the  case,  the  restoration  phase  is 
terminated.  Otherwise,  return  to  (a)  and  continue  the  pro- 
cess until  satisfaction  of  (18-1)  occurs. 

9.2.  Conjugate  Gradient  Phase.  This  phase  involves  a 
single  iteration  and  can  be  summarized  as  follows. 

(a)  Assume  nominal  functions  x(t),  u(t),  tt  which  satisfy 
the  constraints  (2)- (6)  within  the  preselected  accuracy  (18-1). 

(b)  For  the  nominal  functions  and  for  >*=0,  solve  the 
LTP-BVP  (51) -(60)  using  the  method  of  particular  solutions. 

In  this  way,  obtain  the  auxiliary  functions  A*(t),  B*(t),  C* 
and  the  multipliers  X*(t),  p*(t),  a*,  u*  • 

(c)  Set  the  directional  coefficient  > at  the  level  (96) 
for  the  first  conjugate  gradient  phase  and  at  the  level  (97) 
for  any  subsequent  conjugate  gradient  phase.  In  the  latter 
case,  accept  the  directional  coefficient  only  if  Ineq.  (99) 


is  satisfied;  otherwise,  reset  y at  the  level  (96). 

(d)  Compute  the  basic  functions  A(t),  B(t),  C and  the 
multipliers  X (t),  p(t),  o , u using  Eqs.  (101)- (102). 

(e)  Using  the  functions  in  (d) , compute  the  conjugate 
gradient  stepsize  by  a one-dimensional  search  on  the  augmented 
functional  J (a) until  satisfaction  of  Ineq.  (106)  occurs.  Then, 
bisect  the  resulting  stepsize  aQ  (if  necessary) , until  satis- 
faction of  Ineqs.  (105)  occurs. 

(f)  Once  the  conjugate  gradient  stepsize  is  known, 
compute  the  varied  functions  x(t),  u(t),  7T  with  Eqs.  (103). 

9.3.  Conjugate  Gradient-Restoration  Cycle.  After  the 
restoration  phase  is  completed,  verify  whether  Ineq.  (107)  is 
satisfied.  If  this  is  the  case,  start  the  next  cycle  of  the 
sequential  conjugate  gradient-restoration  algorithm.  If  not, 
return  to  the  previous  conjugate  gradient  phase  and  reduce 
the  conjugate  gradient  stepsize  (using  a bisection  process) 
until,  after  restoration,  Ineq.  (107)  is  satisfied. 


10 . Experimental  Conditions 


In  order  to  evaluate  the  theory,  twelve  examples  were 
solved.  The  sequential  conjugate  gradient-restoration 
algorithm  was  programmed  in  FORTRAN  IV,  and  the  numerical 
results  were  obtained  in  double-precision  arithmetic. 

Computations  were  performed  at  Rice  University  using  an 
IBM  370/155  computer.  For  each  example,  the  interval  of  in- 
tegration was  divided  into  100  steps.  The  differential  equa- 
tions were  integrated  using  Hamming's  modified  predictor- 
corrector  method  with  a special  Runge-Kutta  starting  procedure 
(Ref.  10).  The  def ini te * integrals  I,  J,  P,  Q were  computed 
using  a modified  Simpson's  rule.  The  method  of  particular 
solutions  (Refs.  6-8)  was  used  to  solve  the  linear,  two-point 
boundary-value  problems  associated  with  both  the  conjugate 
gradient  phase  and  the  restoration  phase. 

10.1.  Convergence  Conditions.  The  parameters  e^,  e2, 

14 

appearing  in  Ineqs.  (18)  and  (106)  were  set  at  the  levels 

£1  = E-08,  e2  = E-04,  e4  = E-  03.  (108) 

The  tolerance  level  (108-1)  characterizes  the  restoration 
phase;  the  tolerance  levels  (108-1)  and  (108-2) , employed  in 


^■4The  symbol  E + ab  stands  for  lO-3*3. 


combination,  characterize  the  algorithm  as  a whole;  and  the 
tolerance  level  (103-3)  characterizes  the  one-dimensional 
search  for  the  conjugate  gradient  stepsize. 

10.2.  Safeguards . For  the  conjugate  gradient  phase, 
the  parameter  P*  appearing  in  Ineq.  (105-2)  was  set  at  the 
level 


P*  = 10.  (109) 

The  tolerance  level  (109)  limits  the  constraint  violation 
which  is  permissible  during  the  conjugate  gradient  phase. 

Also  for  the  conjugate  gradient  phase,  the  number  of  Hermitian 
search  steps  required  to  satisfy  Ineq.  (106)  was  subject  to 
the  upper  bound 


N < 10.  (110) 

s - 

10.3.  Nonconvergence  Conditions.  The  sequential  conju- 
gate gradient-restoration  algorithm  was  programmed  to  stop 
whenever  satisfaction  of  any  of  the  following  inequalities 
occurred : 


(i) 

N > 50  , 

(111) 

(ii) 

N > 30  , 

(112) 

c 

(iii) 

Nr  > 10  , 

(113) 
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(iv) 

N.  MO  , 
br 

(114) 

(v) 

J (0)  "0, 

3 = 0, 

(115) 

a - 

(vi) 

N.  " 10  , 
bg 

Y = 0 , 

(116) 

(vii) 

N.  ' 10  , 
be 

Y = 0 . 

(117) 

More,  N is  the  total  number  of  iterations,  N is  the  number  of 

c 

cycles,  N is  the  number  of  restorative  iterations  per  cycle, 

N,  is  the  number  of  bisections  of  the  restoration  steps ize 
br 

required  to  satisfy  Ineqs.  (50),  is  the  number  of  bisec- 

tions of  the  conjugate  gradient  stepsize  required  to  satisfy 
Ineqs.  (105),  N^c  is  the  number  of  bisections  of  the  conjugate 
gradient  stepsize  required  to  satisfy  Ineq.  (107),  and  J^(0) 
is  the  slope  of  the  augmented  functional,  at  oi  = 0 . 

Inequalities  (111)  — (112)  apply  to  the  algorithm  as  a 
whole.  Satisfaction  of  (111)  and/or  (112)  is  indicative  of 
extreme  slowness  of  convergence. 

Inequalities  (113)  — (114)  apply  to  the  restoration  phase. 

Satisfaction  of  (113)  is  indicative  of  failure  to  produce  a 
feasible  solution  in  a reasonable  number  of  restorative  itera- 
tions. Satisfaction  of  (114)  is  indicative  of  extreme  small- 
ness of  the  restorative  displacements. 

Inequalities  (115) -(116)  apply  to  the  conjugate  gradient 

! i 

! | 

i ^ 


I 
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phase.  Satisfaction  of  (115)  means  that  the  descent  property 
of  the  conjugate  gradient  phase  does  not  hold,  owing  to 
numerical  inaccuracy.  Satisfaction  of  (116)  is  indicative  of 
extreme  smallness  of  the  conjugate  gradient  displacements. 

Inequality  (117)  applies  to  a complete  conjugate 
gradient-restoration  cycle.  Satisfaction  of  (117)  is  indica- 
tive of  extreme  smallness  of  the  displacements  produced  within 
a complete  conjugate  gradient-restoration  cycle. 

10.4.  Restarting  Conditions.  The  directional  coeffi- 
cient \ of  the  present  conjugate  gradient  phase  was  reset  at 
the  level  > = 0 whenever  satisfaction  of  any  of  the  following 
inequalities  occurred: 


(i) 

J (0)  > 0 , 
a - ' 

Y 

= Q/Q  , 

(118) 

(ii) 

%, ' 10  - 

Y 

= q/6  , 

(119) 

(iii) 

Nbc  > 10  ' 

Y 

= Q/Q  • 

(120) 

Satisfaction  of  (118)  means  that  the  descent  property  of 
the  conjugate  gradient  phase  does  not  hold.  Satisfaction  of 
(119)  is  indicative  of  extreme  smallness  of  the  conjugate 
gradient  displacements.  And  satisfaction  of  (120)  is  indica- 
tive of  extreme  smallness  of  the  displacements  produced  within 
a complete  conjugate  gradient-restoration  cycle. 


i 

\ 


.1 

1 


The  directional  coefficient  y of  the  next  conjugate 
gradient  phase  was  reset  at  the  level  y = 0 whenever  satisfac- 
tion of  any  of  the  following  inequalities  occurred: 


(iv) 

1 ; Nbg  - 10  ' 

y = 0 

or 

Y = Q/Q, 

(121) 

(v) 

1 < Nbc ; 10  • 

o 

ll 

>- 

or 

Y = Q/Q. 

(122) 

Satisfaction  of  (121)  or  (122)  is  indicative  of  large  viola- 
tions of  the  orthogonality  and  conjugacy  conditions,  owing  to 
the  fact  that  the  optimal  conjugate  gradient  stepsize  cannot 
be  employed. 
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11 . Numerical  Examples 

In  this  section,  twelve  numerical  examples  are  des- 
cribed employing  scalar  notation.  In  particular,  the  symbols 
x^(t),  i=l,...,n,  denote  the  components  of  the  state;  the 
symbols  u^(t),  i=l,..,m,  denote  the  components  of  the  control; 
and  the  symbols  tt^,  i=l , . . ,p,  denote  the  components  of  the 
parameter . 

For  all  of  the  examples,  a time  normalization  is  used  in 
order  to  simplify  the  numerical  computations.  Specifically, 
the  actual  time  0 is  replaced  by  the  normalized  time 

t = e/T,  (123) 

which  is  defined  in  such  a way  that  t = 0 at  the  initial  point 
and  t = 1 at  the  final  point.  The  actual  final  time  r,  if  it 
is  free,  is  regarded  as  a component  of  the  vector  parameter  n 
to  be  optimized.  In  this  way,  an  optimal  control  problem  with 
variable  final  time  is  converted  into  an  optimal  control  prob- 
lem with  fixed  final  time. 

Example  11.1.  This  is  a problem  with  (i)  free  initial 
state  and  (ii)  fixed  final  time  t = 1 : 

I = j^(xj  -E  x^  t uj  + u^)dt,  (124) 

*1  = x2  ' *2  = Ul"  X2' 


(125) 
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- U2  + 5t  - 3 = 0 , (126) 

X1  (1)  = 0 . (127) 

The  assumed  nominal  functions  are: 

x1(t)=0,  x2(t)=-l,  u1(t)=0,  u2(t)=0.  (128) 

The  numerical  results  are  given  in  Tables  1-2.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 3 iterations, 
which  include  1 restorative  iteration  and  2 conjugate  gradient 
iterations . 

Example  11.2.  This  is  a problem  with  (i)  a linear  rela- 
tion between  the  components  of  the  initial  state  and  (ii) 
fixed  final  time  t = 1: 


1 = \ ^X1  + x2  + ui/200  + u2)dt, 

(129) 

xL  = x2  , x2  = u1  + u2-x2  , 

(130) 

u1  + 2u2  - 10/(1  + lot)  2 = 0, 

(131) 

Xj^  (0)  + x2  (0)  = -1 , 

(132) 

xx  (1)  + x2  (1)  =1  . 

(133) 

The  assumed  nominal  functions  are: 


mam 
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x1(t)=0,  x2(t)=2t-l,  u1(t)=0,  u2(t)=0.  (134) 

The  numerical  results  are  given  in  Tables  3-4.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 3 iterations, 
which  include  1 restorative  iteration  and  2 conjugate  gradient 
iterations . 

Example  11.3.  This  is  a problem  with  (i)  initial  state 
given  and  (ii)  fixed  final  time  t = 1: 


1 = J (x2 + u^)dt  , (135) 

• 2 

Xj^-X  -u1#  x2  = u2  , (136) 

x^  - u^  - 2x2^2 = 0 , (137) 

xx(0)  = 1,  x2(0)  = ✓(0-1) , (138) 

x1 (1)  = 1 . (139) 


This  problem  is  equivalent  (Ref.  2)  to  that  of  minimizing  (135), 
subject  to  (136-1) , (138-1) , (139)  , and  the  following  state 
inequality  constraint: 

xx  - 0.9  > 0 . (140) 

The  assumed  nominal  functions  are: 

x2  (t)  = /(0. 1)  , u1(t)  = l,  u2(t)=l. 


x^ (t)  = 1, 


(141) 
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The  numerical  results  are  given  in  Tables  5-6.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 14  iterations, 
which  include  9 restorative  iterations  and  5 conjugate  gra- 
dient iterations. 

Example  11.4.  This  is  a problem  with  (i)  initial  state 
partially  given  and  (ii)  fixed  final  time  x = l: 


■£[ 


2x2ul/  (1 


“M 


dt  + *2  (0) 


(142) 

(143) 


u 


1 


(144) 


x1(0)  = 0, 


(145) 


xx (1)  = 1/3  . 


(146) 


This  problem  is  equivalent  (Ref.  2)  to  that  of  minimizing 
(142),  subject  to  (143) , (145),  (146),  and  the  following  con- 
trol inequality  constraint: 

u1  > 0 . (147) 

The  assumed  nominal  functions  are:"^ 


15 


Here,  C=/(l/7).  For  this  value  of  C,  the  constraints 
(143)-(146)  are  satisfied. 
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xx(t)=  (C2/3)(t3  + 3t2  + 3t)  , x2(t)  = C(t  + 1)  , 


(t)  = C,  u2  (t)  = /C. 


(148) 


(149) 


The  numerical  results  are  given  in  Tables  7-8.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 7 iterations, 
which  include  4 restorative  iterations  and  3 conjugate  gra- 
dient iterations. 

Example  11,5.  This  is  a minimum  time  problem  with  (i)  a 
linear  relation  between  the  components  of  the  initial  state 
and  (ii)  free  final  time  t.  After  setting  tt^  = r , the  problem 
is  as  follows: 


I = ii  ± , 


X1  = niui ' 


2 2 


x2  = n x (u^  - xL  - 1/2)  , 


2 2 2 . 
U1  ~ X1  ‘ u2  = °' 


Xj (0) - x2 (0) = 0 , 


xL (1)  = 1,  x2(l)=  -n/4  . 


(150) 


(151) 


(152) 


(153) 


(154) 


This  problem  is  equivalent  (Ref.  2)  to  that  of  minimizing 
(150),  subject  to  (151),  (153),  (154),  and  the  following 
state-derivative  inequality  constraint: 
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x2/n1  + 1/220.  (155) 

The  assumed  nominal  functions  are: 

x1(t)=t,  x2  (t)  = - (7i/4)  t,  (156) 

ux(t)  = 1,  u2(t)  = 1,  tt1=  1.  (157) 

The  numerical  results  are  given  in  Tables  9-10.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 18  iterations, 
which  include  13  restorative  iterations  and  5 conjugate  gra- 
dient iterations. 

Example  11.6.  This  is  a problem  with  (i)  a component  of 
the  initial  state  given,  (ii)  a nonlinear  relation  concerning 
the  remaining  component  of  the  initial  state,  and  (iii)  fixed 
final  time  x = 1 : 

I^j1  |^2x2ux/(l+ux)j  dt  + x“  (0)  , (158) 

x - x*  , x2  = ui  , (159) 

ux  - u2  = 0 , (160) 

x1(0)=0,  j^x2  (0)  - 0.35]  [l  - x2  (0)]  - v*  = 0 , (161) 


xx(l)  = 1/2  . 


(162) 
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This  problem  is  equivalent  (Ref.  2)  to  that  of  minimizing 
(158),  subject  to  (159),  (161-1),  (162),  and  the  following 
inequality  constraints  imposed  on  the  control  and  the  initial 
state : 


u.  > 0, 


(163) 


0 . 35  < *2  (0)  < 1 . 


(164) 


The  assumed  nominal  functions  are: 


xL(t)  = (t3  + 3t2  + 3t)/14,  x2  (t)  = C(t  + 1),  (165) 


u1(t)=C,  u2(t)=/C,  tt1=  /[  (C-  0.35)(1-C)]  . (1661 


The  numerical  results  are  given  in  Tables  11-12.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 10  iterations, 
which  include  7 restorative  iterations  and  3 conjugate  gra- 
dient iterations. 

Example  11.7.  This  is  a problem  with  (i)  a nonlinear 
relation  between  the  components  of  the  initial  state  and  (ii) 
fixed  final  time  x = 1: 


’Here,  C = /(1/7) . 


i. 

c 
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1. 

c 

L 

C 

c 


*• 


L 

L 

0 

c 

l 


9m 

L 

[ 

f 

I 

E 

0 


I = l1  (x2  + u2)dt,  (167) 

J0 

• 2 

xl"xl“ui.  x2  = u2'  (168) 

x2  - - 1 + 2t  - 2x2u2  = 0,  (169) 

x^O)  - x2(0)  = 0.8,  (170) 

x1(l)  = 1 . (171) 


This  problem  is  equivalent  (Ref.  2)  to  that  of  minimizing 
(167) , subject  to  (168-1) , (171) , and  the  following  state 
inequality  constraint: 

x1-0.8-t  + t2>0  . (172) 

The  assumed  nominal  functions  are: 

x1(t)=l,  x (t)  = /(0 . 2)  , u1(t)=l,  u2(t)=0.  (173) 

The  numerical  results  are  given  in  Tables  13-14.  Convergence 
to  the  desired  stopping  condition  occurs  in  N=22  iterations, 
which  include  13  restorative  iterations  and  9 conjugate  gra- 
dient iterations. 

Example  11.8.  This  is  a problem  with  (i)  a component  of 


the  initial  state  given,  (ii)  a linear  relation  between  the 
remaining  components  of  the  initial  state,  and  (iii)  fixed 
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final  time  t = tt/2: 


(174) 


• 2 • 

*l=Tul'  x2=T(2-4xp,  x3  = xx4,  x4  = tu2,  (175) 

2 

4Xiui  - x3u2  - x4  = 0,  (176) 

xx(0)=0,  2x3(0)  + x4(0)  = 1,  (177) 

x1(l)  = 1,  x2(l)  = 0.  (178) 


The  assumed  nominal  functions  are: 

X;L(t)=t,  x2(t)  = 4t(l  - t)  , x3(t)=l-2t,  (179) 

x4 ( t)  = -1 , u1(t)=l/x,  u2(t)=0.  (180) 

The  numerical  results  are  given  in  Tables  15-16.  Convergence 
to  the  desired  stopping  condition  occurs  in  N=23  iterations, 
which  include  13  restorative  iterations  and  10  conjugate  gra- 
dient iterations. 

Example  11.9.  This  is  a problem  with  (i)  a component  of 
the  initial  state  given,  (ii)  a nonlinear  relation  between 
the  remaining  components  of  the  initial  state,  and  (iii)  fixed 
final  time  t = tt/2: 
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I = ! T (uj  - x^  + Tt)dt,  (181) 


2 

*1  = TU1,  x2  = T (2  - 4x1) , 

x3  = tx4 ' 

x4 = iu2. 

(182) 

4xlui  - X3U2  ‘ 

- X4  = 0 , 

(183) 

XL  (0)  =0,  x2(0)  +x3 

(0)  [x3(0)  + 2x4 

(0)]  = -1, 

(184) 

x1(l)  = 1, 

x2(l)  = 0. 

(185) 

The  assumed  nominal  functions 

are : 

x1(t)=t,  x2  (t)  = 4t  (1  - t). 

x^ (t) =l-2t , 

x4  (t)=  - 1, 

(186) 

n1 (t) = 1/ T , 

u 2 (t) = 0 . 

(187) 

The  numerical  results  are  given  in  Tables  17-18.  Convergence 

to  the  desired  stopping  condition  occurs  in  N = 20  iterations, 

which  include  12  restorative  iterations  and  8 conjugate  gra- 
dient iterations. 

Example  11.10.  This  is  a problem  with  (i)  a component 

of  the  initial  state  given,  (ii)  two  nonlinear  relations  between 
the  remaining  components  of  the  initial  state,  and  (iii) 
fixed  final  time  t = tt/2: 
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(188) 


i1  = Tu1,  x2  = T(2-4xp,  x3  = tx4,  x4  = tu2,  (189) 


4x1u1  " x3u2  - x4  = 0, 


(190) 


x1(0)=0,  x2(0)  + x3  (0)  = 1,  x3  (0)x4  (0)  = -1,  (191) 

x3  (1)  =1,  x2  (1)  = 0.  (192) 

This  problem  is  equivalent  (Ref.  2)  to  that  of  minimizing 
(188),  subject  to  ( 189-1 T , (189-2),  (191-1),  (192),  and  the 

following  state  inequality  constraint: 


1 - x2  > 0 . 

(193) 

The  assumed 

nominal 

functions 

are : 

x3 (t) = t , 

x2 (t) = 

4t  ( 1-t), 

x 3 ( t ) * ' 1 -2t,  x4  (t)=  -1 , 

(194) 

u ( t ) 

= 1/T, 

u2 ( t) = 0 . 

(195) 

The  numerical  results  are  given  in  Tables  19-20.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 18  iterations, 
which  include  11  restorative  iterations  and  7 conjugate  gra- 


J 


- J 


0 


.1 

:j 


J 

J 

] 


dient  iterations. 

Example  11.11.  This  is  a minimum  time  problem  with  (i) 
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a nonlinear  relation  between  the  components  of  the  initial 
state  and  (ii)  free  final  time  x.  After  setting  ti^  = t,  the 
problem  is  as  follows: 


rH 

t= 

II 

M 

(196) 

x2  = XT  1 (u^  - x^  - 1/2)  , 

(197) 

2 2 2 

U1  ~ X1  ~ u2  = 0 ' 

(198) 

xL(0)  + x2(0)  = 0, 

(199) 

% 

x1  (1)  x2(l)  =-tt/4  . 

(200) 

This  problem  is  equivalent  (Ref.  2)  to  that  of  minimizing 
(196) , subject  to  (197) , (199)  , (200)  , and  the  following 
state-derivative  inequality  constraint: 

x2/n1 + 1/2  > 0 . (201) 

The  assumed  nominal  functions  are: 

xx(t)=t,  x2(t)=-(TT/4)t,  u1(t)=l,  u2(t)=l,  n1=l.  (202) 

The  numerical  results  are  given  in  Tables  21-22.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 14  iterations, 
which  include  10  restorative  iterations  and  4 conjugate  gra- 


dient iterations. 
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Example  11.12.  This  is  a problem  with  (i)  a linear 
relation  between  the  first  two  components  of  the  initial 
state,  (ii)  a nonlinear  relation  between  the  remaining  compo- 
nents of  the  initial  state,  and  (ii)  fixed  final  time  t = 1: 

1 = 1 u^dt , (203) 

*0 

x1  = x2  , x2  = u^,  x3  = x4 , x4  = u2  , (204) 

+ 2x3u2  + 2x4  = 0,  (205) 

xx  (0)  + x2  (0)  = 1,  x3  (0)  + 2x3  (0)x4  (0)= -0.85,  (206) 

x^l)  = 0,  x2(l)  = -1.  (207) 

The  assumed  nominal  functions  are: 

x1(t)=0,  x2(t)=l-2t,  x3  (t)  = (l-2t) /(0.15),  (208) 

x4(t)=  (2t-l)/2/(0.15),  ux(t)=l,  u2(t)=  0 . (209) 

The  numerical  results  are  given  in  Tables  23-24.  Convergence 
to  the  desired  stopping  condition  occurs  in  N = 20  iterations, 
which  include  12  restorative  iterations  and  8 conjugate  gra- 
dient iterations. 


■— 


Nc  Ng  Nr  N Y P Q I 


0000  0.43E+01 

r-  1011  0.15E-29  0.17E+00  1.29984 

2102  y=0  0.90E-30  0.49E-02  1.18165 

3103  y^O  0.12E-29  0.31E-06  1.18040 

J1  

L 


r 


Table  2.  Converged  solution.  Example  11.1. 


t 

X1 

X2 

U1 

U2 

0.0 

-0.0069 

-0.3314 

1.5004 

-1.4995 

0.1 

-0.0316 

-0.1699 

1.2368 

-1.2631 

0.2 

-0.0422 

-0.0485 

0.9798 

-1.0201 

0.3 

-0.0425 

0.0370 

0.7279 

-0.7720 

0.4 

-0.0359 

0.0908 

0.4799 

-0.5200 

0.5 

-0.0253 

0.1159 

0.2343 

-0.2656 

0.6 

-0.0135 

0.1153 

-0.0100 

-0.0100 

0.7 

-0.0030 

0.0916 

-0.2544 

0.2455 

0.8 

0.0040 

0.0468 

-0.5000 

0.4999 

0.9 

0.0056 

-0.0171 

-0.7481 

0.7518 

1.0 

0.0000 

-0.0988 

-1.0000 

1.0000 

I 

I I 
[ 


x = 1.00000 
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Table  3 . Convergence 

history. 

Example  11.2. 

N N N N v 

c g r 1 

P 

Q I 

0 

0 

0 

0 

0.79E+01 

1 

0 

1 

1 

0.37E-27 

0. 35E+00 

1.89493 

2 

1 

0 

2 

Y=0 

0.15E-26 

0.65E-02 

1.48434 

3 

1 

0 

3 

0.24E-26 

0.50E-04 

1.48042 

Table  4 . 

Converged  solution,  Example  11.2. 

t 

X1 

x2  u2 

0.0 

0.1127 

-1.1127 

12.0825 

-1.0412 

0.1 

0.0416 

-0.4340 

4.6589 

-1.0794 

0.2 

0.0146 

-0.1335 

3.2253 

-1.0571 

0.3 

0.0114 

0.0578 

2.7069 

-1.0409 

0.4 

0.0245 

0.1984 

2.4732 

-1.0366 

0.5 

0.0501 

0.3101 

2.3668 

-1.0445 

0.6 

0.0859 

0.4035 

2.3329 

-1.0644 

0.7 

0.1304 

0.4847 

2.3482 

-1.0959 

0.8 

0.1826 

0.5581 

2.4013 

-1.1389 

0.9 

0.2419 

0.6265 

2.4864 

-1.1932 

1.0 

0.3078 

0.6921 

2.5999 

-1.2586 

t = 1.00000 


y 

i. 
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Table  5.  Convergence  history.  Example  11.3. 


Nc 

N 

g 

N 

r 

N 

Y 

P 

Q 

I 

0 

0 

0 

0 

0. 14E+01 

1 

0 

3 

3 

0. 52E-09 

0.35E+00 

1.83569 

2 

l 

2 

6 

y=0 

0.15E-16 

0.14E-01 

1.66599 

3 

l 

1 

8 

y^O 

0.12E-08 

0.27E-03 

1.65745 

4 

l 

1 

10 

Y^O 

0.12E-12 

0.18E-03 

1.65704 

5 

l 

1 

12 

Y^O 

0.51E-13 

0.16E-03 

1.65677 

6 

l 

1 

14 

Y^O 

0.34E-12 

0.66E-04 

1.65641 

Table  6.  Converged  solution , Example  11.3. 

t 

xx  x2  u2 

0.0 

1.0000 

0.3162 

1.7514 

-1.1881 

0.1 

0.9417 

0.2043 

1.3437 

-1.1180 

0.2 

0.9082 

0.0907 

1.0166 

-1.0558 

0.3 

0.9002 

0.0151 

0.8232 

-0.4215 

0.4 

0.9000 

-0.0027 

0.8099 

-0.0180 

0.5 

0.9000 

-0.0013 

0.8100 

0.0096 

0.6 

0.9000 

-0.0010 

0.8100 

0.0303 

0.7 

0.9003 

0.0175 

0.7956 

0.4250 

0.8 

0.9087 

0.0933 

0.6292 

1.0524 

0.9 

0.9420 

0.2051 

0.4398 

1.0912 

1.0 

1.0000 

0.3162 

0.2376 

1.2053 

t = 1.00000 

[ 
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Table  7.  Convergence 

history. 

Example  11.4. 

N N N Y 

g r 

P 

Q I 

0 

0 

0 

0.47E-32 

0.49E-01 

0.19642 

1 

2 

3 

Y = 0 

0.10E-13 

0.11E-02 

0.16924 

1 

1 

5 

y/o 

0.12E-09 

0.10E-03 

0.16816 

1 

1 

7 

y?o 

0.31E-13 

0.11E-04 

0.168U9 

.1 


Table  8.  Converged  solution.  Example  11.4. 


t 

X1 

x2 

U1 

U2 

0.0 

0.0000 

0.2908 

0.9570 

0.9782 

0.1 

0.0112 

0.3753 

0.7474 

0.8645 

0.2 

0.0281 

0.4428 

0.6127 

0.7827 

0.3 

0.0504 

0.4995 

0.5266 

0.7257 

0.4 

0.0779 

0.5488 

0.4603 

0.6784 

0.5 

0.1105 

0.5917 

0.3978 

0.6307 

0.6 

0.1478 

0.6284 

0.3372 

0.5807 

0.7 

0.1894 

0.6593 

0.2817 

0.5308 

0.8 

0.2346 

0.6848 

0.2269 

0.4764 

0.9 

0.2829 

0.7039 

0.1491 

0.3862 

1.0 

0.3333 

0.7138 

0.0516 

0.2272 

r = 1.00000 
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Table  9.  Convergence  history.  Example  11.5. 


Nc 

N 

g 

N 

r 

N 

y 

P 

Q 

I 

0 

0 

0 

0 

0.11E+01 

1 

0 

5 

5 

0.13E-14 

0.21E-01 

1.83056 

2 

i 

2 

8 

>=o 

0.52E-14 

0.37E-02 

1.82421 

3 

i 

2 

11 

Y^O 

0.18E-15 

0.20E-02 

1.82301 

4 

i 

2 

14 

Y/0 

0.72E-15 

0.11E-02 

1.82249 

5 

l 

1 

16 

y/0 

0.11E-08 

0.27E-03 

1.82230 

6 

i 

1 

18 

Y^O 

0. 81E-11 

0.71E-04 

1.82226 

Table  10. 

Converged  solution 

l.  Example 

11.5. 

t 

X1 

x2 

U1 

U2 

0.0 

0.0034 

0.0034 

0.4978 

0.4978 

0.1 

0.0936 

-0,0436 

0,4887 

0.4797 

0.2 

0.1805 

-0.0967 

0.4631 

0.4264 

0.3 

0.2614 

-0.1609 

0.4225 

0.3319 

0.4 

0.3338 

-0.2396 

0.3703 

0.1601 

0.5 

0.4020 

-0.3298 

0.4025 

-0.0193 

0.6 

0.4824 

-0.4209 

0.4824 

0.0013 

0.7 

0.5788 

-0.5120 

0.5788 

-0.0015 

0.8 

0.6945 

-0.6031 

0.6945 

0.0027 

0.9 

0.8334 

-0.6942 

0.8334 

-0.0029 

1.0 

1.0000 

-0.7853 

1.0000 

-0.0016 

T = tt1  = 1.82226 

' 
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Table 

11. 

Convergence  history, 

Example  11 

.6. 

N 

c 

N 

g 

N 

r 

N 

y 

p 

Q 

I 

1" 

0 

0 

0 

0 

0. 31E-01 

1 

0 

3 

3 

0. 28E-09 

0.65E-01 

0.31197 

2 

l 

2 

6 

7=0 

0.14E-15 

0.17E-02 

0.28910 

3 

l 

1 

8 

y?o 

0.44E-10 

0.27E-03 

0.28781 

4 

l 

1 

10 

y/o 

0.52E-12 

0.19E-04 

0.28762 

Table  12.  Converged  solution,  Example  11.6. 


0.0 

0.0000 

0.4133 

0.9597 

0.9796 

0.1 

0.0210 

0.4990 

0.7660 

0.8752 

0.2 

0.0496 

0.5686 

0.6354 

0.7971 

0.3 

0.0855 

0.6275 

0.5471 

0.7396 

0.4 

0.1283 

0.6787 

0.4789 

0.6920 

0.5 

0.1775 

0.7235 

0.4174 

0.6460 

0.6 

0.2328 

0.7623 

0.3576 

0.5980 

0.7 

0.2935 

0.7950 

0.2974 

0.5454 

0.8 

0.3590 

0.8216 

0.2315 

0.4812 

0.9 

0.4282 

0.8409 

0.1525 

0.3906 

1.0 

0.5000 

0.8518 

0.0677 

0.2602 

1 = 1.00000  , ir1  = 0.19278 
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Table  13.  Convergence  history,  Example  11.7. 


N 

c 

N 

g 

N 

r 

N 

7 

P 

Q 

I 

0 

0 

0 

0 

0. 33E+00 

1 

0 

3 

3 

0.52E-14 

0.40E+01 

2.57561 

2 

1 

2 

6 

7*0 

0.32E-11 

0.11E-01 

1.82401 

3 

1 

1 

8 

77*0 

0. 25E-09 

0.20E-01 

1.81492 

4 

1 

1 

10 

75*0 

0.18E-09 

0.18E-02 

1.80922 

5 

1 

1 

12 

77*0 

0.11E-11 

0.22E-02 

1.80803 

6 

1 

1 

14 

77*0 

0.31E-12 

0.40E-03 

1.80736 

7 

1 

1 

16 

77*0 

0.37E-13 

0.74E-03 

1.80704 

8 

1 

1 

18 

7 ^0 

0.35E-13 

0.15E-03 

1.80678 

9 

1 

1 

20 

77*0 

0.15E-14 

0. 32E-03 

1.80667 

10 

1 

1 

22 

77*0 

0.14E-13 

0.98E-04 

1.80650 

Table  14.  Converged  solution,  Example  11.7. 

t 

X1  x2  U1  u2 

0.0 

0.8139 

0.1180 

0.0102 

-1.4724 

0.1 

0.8901 

0.0124 

0.0056 

-0.5342 

0.2 

0.9600 

-0.0063 

0.3219 

0.0235 

0.3 

1.0100 

0.0000 

0.6201 

0.0633 

0.4 

1.0400 

0.0033 

0.8816 

0.0012 

0.5 

1.0500 

-0.0001 

1.1024 

-0.0725 

0.6 

1.0400 

-0.0085 

1.2808 

-0.0525 

0.7 

1.0101 

0.0108 

1.4068 

0.6212 

0.8 

0.9779 

0.1340 

1.1179 

1.6351 

0.9 

0.9747 

0.2910 

0.8558 

1.5361 

1.0 

1.0000 

0.4472 

0.6327 

1.5286 

t = 1.00000 
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Table  15.  Convergence  history.  Example  11.8 


N 

c 

N 

g 

N 

r 

N 

Y P 

Q 

I 

0 

0 

0. 23E+01 

1 

3 

3 

0.49E-09 

0. 52E-01 

1.27196 

2 

l 

1 

5 

Y=0  0.86E-09 

0.27E-01 

1.25174 

3 

l 

1 

7 

>70  0.83E-09 

0.10E-01 

1.24337 

4 

l 

1 

9 

Y5*0  0.54E-10 

0.45E-02 

1.23945 

5 

l 

1 

11 

Y^O  0.85E-11 

0.32E-02 

1.23734 

6 

l 

1 

13 

Yt*0  0.40E-11 

0. 22E-02 

1.23589 

7 

l 

1 

15 

YYO  0.13E-11 

0.13E-02 

1.23492 

8 

l 

1 

17 

Y^o  0.22E-12 

0.79E-03 

1.23435 

9 

l 

1 

19 

Y^0  0.31E-13 

0.43E-03 

1.23403 

10 

l 

1 

21 

YY0  0.33E-14 

0.21E-03 

1.23385 

11 

l 

1 

23 

Y^O  0.24E-15 

0.96E-04 

1.23377 

Table 

16. 

Converged  solution,  Example  11.8 

• 

B 1 

> 

C1 

> 

c2  x3  x4 

U1 

U2 

0.0 

0.0000 

-0. 

0122  0.9981  -0.9963 

0.9955 

-0.9944 

0.1 

0.1557 

0. 

2968  0.8300  -1.1396 

0.9715 

-0.8354 

0.2 

0.3076 

0. 

5758  0.6416  -1.2548 

0.9557 

-0.6208 

0.3 

0.4516 

0. 

7978  0.4376  -1.3403 

0.8747 

-0.4941 

0.4 

0.5850 

0. 

9417  0.2221  -1.3953 

0.8130 

-0.2003 

0.5 

0.7044 

0. 

9927  0.0014  -1.4094 

0.7050 

0.0059 

0.6 

0.8062 

0. 

9464  -0.2191  -1.3944 

0.5899 

0.1908 

0.7 

0.8890 

0. 

8070  -0.4350  -1.3478 

0.4611 

0.4060 

0.8 

0.9502 

0. 

5873  -0.6407  -1.2664 

0.3157 

0.6301 

0.9 

0.9876 

0. 

3090  -0.8310  -1.1507 

0.1589 

0.8379 

1.0 

1.0000 

0. 

0000  -1.0007  -1.0047 

-0.0015 

1.0148 

r = tt/2  = 1.57079 


Table  17.  Convergence  history.  Example  11.9 


N 

c 

N 

g 

Nr 

N 

Y 

p 

Q 

I 

0 

0 

0 

0 

0.23E+01 

1 

0 

4 

4 

0.16E-13 

0.22E-01 

1.24364 

2 

l 

1 

6 

Y = 0 

0.24E-10 

0.55E-02 

1.23682 

3 

l 

1 

8 

y/o 

0.20E-11 

0.13E-02 

1.23518 

4 

l 

1 

10 

y/o 

0.28E-13 

0.54E-03 

1.23465 

5 

l 

1 

12 

yyo 

0.20E-14 

0.37E-03 

1.23438 

6 

l 

1 

14 

y/o 

0.55E-15 

0. 28E-03 

1.23419 

7 

l 

1 

16 

y/0 

0.61E-15 

0.20E-03 

1.23403 

8 

l 

1 

18 

Y7*0 

0.93E-15 

0.12E-03 

1.23394 

9 

l 

1 

20 

y/o 

0.44E-15 

0.72E-04 

1.23388 

Table  18. 

Converged  solution,  Example  11.9. 

t 

X1  x2 

x3  x4  ux 

u2 

0.0 

0.0000 

0.0422 

1.0522 

-1.0213 

1.0068 

-0.9914 

0.1 

0.1578 

0.3511 

0.8802 

-1.1643 

1.0051 

-0.8188 

0.2 

0.3122 

0.6292 

0.6879 

-1.2786 

0.9674 

-0.6198 

0.3 

0.4579 

0.8486 

0.4800 

-1.3652 

0.9013 

-0.4433 

0.4 

0.5940 

0.9872 

0.2613 

-1.4117 

0.8199 

-0.1709 

0.5 

0.7142 

1.C306 

0.0382 

-1.4239 

0.7099 

0.0156 

0.6 

0.8166 

0.9746 

-0.1843 

-1.4045 

0.5904 

0.2399 

0.7 

0.8985 

0.8243 

-0.4010 

-1.3469 

0.4496 

0.4938 

0.8 

0.9570 

0.5949 

-0.6054 

-1.2504 

0.2931 

0.7286 

0.9 

0.9906 

0.3106 

-0.7921 

-1.1208 

0.1347 

0.9115 

1.0 

1.0000 

0.0000 

-0.9563 

-0.9674 

-0.0124 

1.0307 

t = tt/2  = 1.57079 
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Table  19.  Convergence  history.  Example  11.10. 


N N N N 
c g r 


0 

0 

0 

0. 23E+01 

0 

4 

4 

0.21E-15 

0.24E-01 

1.24401 

1 

1 

6 

y=0 

0.29E-10 

0.60E-02 

1.23670 

1 

1 

8 

y/0 

0.25E-11 

0.15E-02 

1.23495 

1 

1 

10 

Y/0 

0.35E-13 

0.60E-03 

1.23438 

1 

1 

12 

Y^0 

0.26E-14 

0.40E-03 

1.23408 

1 

1 

14 

Y^0 

0.54E-15 

0.24E-03 

1.23389 

1 

1 

16 

Y^0 

0.55E-16 

0.12E-03 

1.23378 

1 

1 

18 

Y^0 

0.26E-17 

0.48E-04 

1.23374 

Table  20.  Converged  solution.  Example  11.10. 


t 

X1 

X2 

x3 

x4 

U1 

u2 

0.0 

0.0000 

-0.0008 

1.0004 

-0.9995 

1.0016 

-0.9987 

0.1 

0.1563 

0.3081 

0.8317 

-1.1434 

0.9848 

-0.8312 

0.2 

0.3087 

0.5870 

0.6426 

-1.2594 

0.9432 

-0.6553 

0.3 

0.4534 

0.8085 

0.4374 

-1.3456 

0.9002 

-0.4064 

0.4 

0.5878 

0.9506 

0.2221 

-1.3909 

0.8044 

-0.1941 

0.5 

0.7061 

0.9999 

0.0018 

-1.4090 

0.7028 

-0.0310 

0.6 

0.8081 

0.9520 

-0.2189 

-1.3979 

0.5922 

0.1808 

0.7 

0.8910 

0.8104 

-0.4353 

-1.3506 

0.4601 

0.4227 

0.8 

0.9517 

0.5887 

-0.6413 

-1.2655 

0.3100 

0.6568 

0.9 

0.9880 

0.3092 

-0.8311 

-1.1462 

0.1527 

0.8544 

1.0 

1.0000 

0.0000 

-1.0000 

-0.9999 

0.0007 

0.9970 

t = tt/2  = 1.57079 
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Table  21  . Convergence  history.  Example  11.11. 


Nc 

N 

g 

Nr 

N 

Y 

P 

Q 

I 

0 

0 

0 

0 

0.11E+01 

1 

0 

4 

4 

0.28E-13 

0.10E+00 

1.61985 

2 

l 

2 

7 

y=0 

0.39E-12 

0.14E-01 

1.58737 

3 

l 

2 

10 

Y*0 

0.99E-16 

0.14E-02 

1.58313 

4 

l 

1 

12 

Y^O 

0.19E-08 

0.26E-03 

1.58236 

5 

l 

1 

14 

Y/0 

0.10E-09 

0.64E-04 

1.58221 

Table  22.  Converged  solution.  Example  11.11. 

t 

X1  x2  U1  u2 

0.0 

-0.0928 

-0.3047 

0.8187 

0.8134 

0.1 

0.0372 

-0.2771 

0.8228 

0.8219 

0.2 

0.1664 

-0.2526 

0.8066 

0.7893 

0.3 

0.2914 

-0.2414 

0.7704 

0.7131 

0.4 

0.4092 

-0.2526 

0.7144 

0.5856 

0.5 

0.5165 

-0.2930 

0.6402 

0.3781 

0.6 

0.6141 

-0.3627 

0.6184 

0.0732 

0.7 

0.7194 

-0.4417 

0.7195 

-0.0104 

0.8 

0.8428 

-0.5208 

0.8428 

0.0051 

0.9 

0.9872 

-0.5999 

0.9873 

-0.0034 

1.0 

1.1565 

-0.6790 

1.1565 

-0.0087 

t = = 1.58221 

eii& 
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Table  23.  Convergence  history,  Example  11.12. 


Nc  Ng  Nr  N y P Q I 


0 

0 

0 

0 

0.22E+02 

1 

0 

4 

4 

0.17E-08 

0.46E+00 

3.32348 

2 

1 

2 

7 

Y=0 

0.67E-13 

0.55E-01 

3.05708 

3 

1 

2 

10 

Y*0 

0.96E-16 

0.15E-01 

3.01500 

4 

1 

1 

12 

Y=0 

0.11E-16 

0.25E-02 

3.00374 

5 

1 

1 

14 

Y=0 

0.20E-17 

0.13E-02 

3.00161 

6 

1 

1 

16 

Y=0 

0.17E-16 

0.51E-03 

3.00073 

7 

1 

1 

18 

Y=0 

0.31E-19 

0.27E-03 

3.00034 

8 

1 

0 

19 

Y=0 

0.61E-08 

0.11E-03 

3.00016 

9 

1 

0 

20 

Y=0 

0.63E-08 

0.60E-04 

3.00007 

Table  24.  Converged  solution,  Example  11.12. 


t 

X1 

X2 

X3 

X4 

U1 

U2 

0.0 

0.5018 

0.4981 

0.5645 

-1.0351 

-3.0237 

0.7798 

0.1 

0.5370 

0.2121 

0.4651 

-0.9487 

-2.7051 

0.9725 

0.2 

0.5452 

-0.0432 

0.3756 

-0.8348 

-2.4024 

1.3420 

0.3 

0.5294 

-0.2684 

0.2998 

-0.6704 

-2.1021 

2.0061 

0.4 

0.4925 

-0.4637 

0.2443 

-0.4231 

-1.8023 

2.9551 

0.5 

0.4377 

-0.6282 

0.2179 

-0.0969 

-1.4856 

3.3646 

0.6 

0.3679 

-0.7613 

0.2240 

0.2027 

-1.1815 

2.4533 

0.7 

0.2864 

-0.8650 

0.2543 

0.3824 

-0.8926 

1.1796 

0.8 

0.1959 

-0.9396 

0.2968 

0.4533 

-0.5988 

0.3164 

0.9 

0.0994 

-0.9847 

0.3428 

0.4582 

-0.3017 

-0.1722 

1.0 

0.0000 

-1.0000 

0.3872 

0.4253 

-0.0033 

-0.4628 

t = 1.00000 
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12.  Discussion 


The  examples  presented  in  Section  11  were  solved  with 


both  the  sequential  conjugate  gradient-restoration  algorithm 


(SCC.RA)  of  this  report  and  the  sequential  ordinary  gradient- 


restoration  algorithm  (SOGRA)  of  Ref.  5.  This  was  done  in 


order  to  gain  perspective  on  the  relative  merit  of  SCGRA 


vis-a-vis  SOGRA. 


The  comparative  results  are  presented  in  Tables  25-27, 


where  the  number  of  iterations  N required  to  achieve  different 


tolerance  levels  for  the  error  in  the  optimality  conditions  Q 


is  given  for  a fixed  tolerance  level  of  the  constraint  error 


r v E-08.  Also  shown  in  the  tables  are  the  values  obtained 


for  the  objective  functional  I. 


Cumulative  results  for  the  twelve  examples  investigated 


are  given  in  Table  28.  Here,  the  total  number  of  iterations 


for  convergence  )!N  is  presented  as  a function  of  the  tolerance 


level  chosen  for  the  error  in  the  optimality  conditions  Q,  for 


a fixed  tolerance  level  in  the  constraint  error  P<  E-08.  In 


this  comparative  study,  tolerance  levels  in  the  range  0*-  E-02 


to  0 1 E-04  were  chosen  for  the  error  in  the  optimality  conditions. 


In  Tables  25-27,  the  symbol  LQ  stands  for  a linear-quadratic 
problem,  and  the  symbol  NLQ  stands  for  a nonlinear  and/or 
nonquadratic  problem. 


From  Tables  25-28,  it  appears  that  there  is  no  saving 
in  number  of  iterations  for  Q<E-02.  The  saving  is  4% 
for  Q < E-03  and  11%  for  Q<E-04.  Clearly,  the  relative 
advantage  of  SCGRA  with  respect  to  SOGRA  increases  by  imposing 
a tighter  tolerance  level  on  the  error  in  the  optimality  con- 
ditions . 

It  must  be  noted  that  the  experiments  performed  show 
that  the  computer  time  per  iteration  is  roughly  the  same  for 
SCGRA  and  SOGRA.  Therefore,  the  conclusions  pertaining  to 
savings  in  number  of  iterations  also  apply  to  savings  in  com- 
puter time. 
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Table  25. 

Results  for  P < E-08  and  Q<E-02. 

Example  Type 

SCGRA  SOGRA 

N I N I 

11.1 

LQ 

2 

1.18165 

2 

1.18165 

11.2 

LQ 

2 

1.48434 

2 

1.43434 

11.3 

NLQ 

8 

1.65745 

8 

1.65742 

11.4 

NLQ 

3 

0.16924 

3 

0.16924 

11.5 

NLQ 

8 

1.82421 

8 

1.82421 

11.6 

NLQ 

6 

0.28910 

6 

0.28910 

11.7 

NLQ 

10 

1.80922 

10 

1.81240 

11.8 

NLQ 

9 

1.23945 

7 

1.24498 

11.9 

NLQ 

6 

1.23682 

6 

1.23682 

11.10 

NLQ 

6 

1.23670 

6 

1.23670 

11.11 

NLQ 

10 

1.58313 

9 

1.58381 

11.12 

NLQ 

3.00374 

11 

3.00785 

Table  26. 

Results 

for  P < E-08  and  Q<E-03. 

Example  Type 

SCGRA  SOGRA 

N 

M 

z 

w 

n.i 

LQ 

3 

1.18040 

3 

1.18041 

11.2 

LQ 

3 

1.48042 

3 

1.48045 

11.3 

NLQ 

8 

1.65745 

8 

1.65742 

11.4 

NLQ 

5 

0.16816 

5 

0.16818 

11.5 

NLQ 

16 

1.82230 

12 

1.82276 

11.6 

NLQ 

8 

0.28781 

8 

0.28774 

11.7 

NLQ 

14 

1.80736 

18 

1.80797 

11.8 

NLQ 

17 

1.23435 

19 

1.23613 

11.9 

NLQ 

10 

1.23465 

10 

1.23492 

11.10 

NLQ 

10 

1.23438 

11 

1.23439 

11.11 

NLQ 

12 

1.58236 

13 

1.58242 

11.12 

NLQ 

16 

3.00073 

17 

3.00075 

. 


; 
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Example 

Type 

SCGRA 

SOGRA 

N 

I 

N I 

11.1 

LQ 

3 

1.18040 

3 

1.18041 

11.2 

LQ 

3 

1.48042 

3 

1.48045 

11.3 

NLQ 

14 

1.65641 

12 

1.65678 

11.4 

NLQ 

7 

0.16809 

7 

0.16810 

11.5 

NLQ 

18 

1.32226 

20 

1.82224 

11.6 

NLQ 

10 

0.28762 

10 

0.28764 

11.7 

NLQ 

22 

1.80650 

26 

1.80654 

11.8 

NLQ 

23 

1.23377 

35 

1.23400 

11.9 

NLQ 

20 

1.23388 

20 

1.23411 

11.10 

NLQ 

18 

1.23374 

20 

1.23384 

11.11 

NLQ 

14 

1.58221 

16 

1.58217 

11.12 

NLQ 

20 

3.00007 

21 

3.00007 

Table  28. 

Cumulative  number  of 
for  convergence,  P ^ E 

iterations 

-08. 

Q 

SCGRA 

SOGRA 

EN 

EN 

Q < E-02 

82 

78 

Q < E-03 

122 

127 

Q < E-04 

172 

193 

L 


i. 


T 

<* 

1 

! 

T 

*. 

r 

i. 

L 

I 

I 

I 

I 

1 

I 

I 
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13 . Conclusions 

In  this  report,  a new  member  of  the  family  of  sequential 
gradient-restoration  algorithms  for  the  solution  of  optimal 
control  problems  is  presented.  This  is  an  algorithm  of  the 
conjugate  gradient  type  and  solves  the  problem  represented  by 
Eqs.  ( 1 ) - (6 ) : Minimize  a functional  subject  to  differential 
constraints,  nondifferential  constraints,  and  general  boundary 
conditions . 

The  algorithm  presented  here  differs  from  those  of  Refs. 
3-4,  in  that  it  is  not  required  that  the  state  vector  be  given 
at  the  initial  point.  Instead,  the  initial  conditions  can  be 
absolutely  general.  In  analogy  with  Refs.  3-4,  the  present 
algorithm  is  capable  of  handling  general  final  conditions; 
therefore,  it  is  suitable  for  the  solution  of  optimal  control 
problems  with  general  boundary  conditions. 

The  importance  of  the  present  algorithm  lies  in  that  many 
optimal  control  problems  either  arise  naturally  in  the  present 
format  or  can  be  brought  to  such  a format  by  means  of  suitable 
transformations  (see  Ref.  2).  Therefore,  a great  variety  of 
optimal  control  problems  can  be  handled,  as  it  is  shown  by  the 
numerical  examples  presented. 
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* 


a 
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Twelve  numerical  examples  are  presented  to  illustrate  the 
performance  of  the  algorithm.  The  numerical  results  show  the 
feasibility  as  well  as  the  convergence  characteristics  of  the 
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algorithm.  A comparative  analysis  of  the  sequential  conjugate 
gradient-restoration  algorithm  and  the  sequential  ordinary 
gradient-restoration  algorithm  shows  that  the  relative  advan- 
tage of  SCGRA  with  respect  to  SOGRA  increases  by  imposing  a 
tighter  tolerance  level  on  the  error  in  the  optimality  condi- 
tions . 

In  summary,  the  new  member  of  the  family  of  sequential 
gradient-restoration  algorithms  described  here  has  the  follow- 
ing properties:  (i)  it  retains  the  robustness,  reliability, 
and  convergence  characteristics  of  the  algorithms  discussed  in 
Refs.  3-4;  (ii)  it  is  able  to  handle  all  of  the  optimal  con- 
trol problems  treated  in  Refs.  3-4;  and  (iii)  it  has  the  ad- 
ditional capability  of  handling  optimal  control  problems  with 
general  boundary  conditions. 
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